function [slm,xp,yp] = slmengine(x,y,varargin)
% slmengine: estimates a spline function from data plus a fit prescription
% usage 1: slm = slmengine(x,y);
% usage 2: slm = slmengine(x,y,prescription);
% usage 3: slm = slmengine(x,y,prop1,val1,prop2,val2,...);
% usage 4: slm = slmengine(x,y,prescription,prop1,val1,prop2,val2,...);
% usage 5: [slm,xp,yp] = slmengine(x,y,prescription,prop1,val1,prop2,val2,...);
%
% Note: slmengine is the command driven tool for fitting a model
%   to data. It uses either a prescription structure (as supplied by
%   slmset) or sets of property/value pairs. Those pairs are defined
%   in the help to slmset. slmengine is also used by slmfit, the gui
%   tool for fitting a curve to data.
%
% Note: The optimization toolbox (lsqlin) is generally required for
% most fits. Some models (those where the break points are also
% estimated) require fmincon.
%
% Note: Some prescriptive parameters or combinations will be
%   inappropriate for the lower degree models.
%
% arguments: (input)
%  x,y     - vectors of data used to fit the model. Any
%            NaN or inf elements (in x or y) will be excluded
%            from the fit.
%
%  prescription - structure generated by slmset to control the
%            fitting process.
%
%            Alternatively, one can simply supply property/value
%            pairs directly to slmengine.
%
% Arguments: (output)
%  slm     - model as a shape-language-model structure,
%            normally constructed by slmfit or slmengine. slm
%            may also be returned in other forms, either a pp
%            structure tht ppval can use, or as a simple array
%            of coefficients.
%
%            Upon return, slm will contain a structure field 'stats'
%            that holds fields which describe the quality of fit:
%
%            slm.stats.TotalDoF = Total degrees of freedom in the model
%
%            slm.stats.NetDoF = Net degrees of freedom in the spline model.
%               Thus, NetDoF reflects any equality constraints in the model
%
%            slm.stats.R2 = Traditional R-squared coefficient
%
%            slm.stats.R2Adj = Adjusted R^2, accounting for changing
%               degrees of freedom. For a definition, look here:
%
%               http://en.wikipedia.org/wiki/Coefficient_of_determination
%
%            slm.stats.RMSE = Root-Mean-Squared-Error
%
%            slm.stats.ErrorRange = 1x2 vector that contains the minimmum
%               and maximum errors, defined as (Yhat - Y)
%
%            slm.stats.Quartiles = 1x2 vector that contains the 25% and 75%
%               error quartiles. This gives a trimmed measure of the
%               error magnitudes, less any outliers.
%
%  xp, yp  - predicted points along the fitted curve. The number of these
%            points is defined by the predictions property. (see slmset)
%            xp will be a list of equally spaced points. 
%
% Example:
%  (See the SLM_tutorial demo for other examples.)
%  x = rand(1,50);
%  y = sin(x*2*pi) + randn(size(x))/10;
% 
%  slm = slmengine(x,y,'knots',0:.2:1,'plot','on','concavedown','on','minvalue',0)
%  slm = 
%            form: 'slm'
%          degree: 3
%           knots: [6x1 double]
%            coef: [6x2 double]
%    prescription: [1x1 struct]
%               x: [50x1 double]
%               y: [50x1 double]
%
% Returned in a form that ppval or fnval can use:
%
%  slm = slmengine(x,y,'knots',0:.2:1,'plot','on','concavedown','on','minvalue',0,'res','pp')
% slm = 
%            form: 'pp'
%          breaks: [0 0.2 0.4 0.6 0.8 1]
%           coefs: [5x4 double]
%          pieces: 5
%           order: 4
%             dim: 1
%    prescription: [1x1 struct]
%
%
% See also: spap2, ppval, fnval
%
% Author: John D'Errico
% e-mail: woodchips@rochester.rr.com
% Release: 1.0
% Release date: 1/2/08


% which form does the shape prescription take? Or do we
% just use our defaults for the prescription?
if nargin<2
  error('SLMENGINE:improperdata','Must supply two vectors to fit a model: x & y')
elseif nargin==2
  prescription = slmset;
elseif (nargin>=3)
  % prescription supplied directly or as property/value pairs,
  % or as a prescription, modified by a few set of pairs.
  % slmset will resolve all of these cases.
  prescription = slmset(varargin{:});
end

% check the data for size, turning it into column vectors
x = x(:);
y = y(:);
n = length(x);
if n~=length(y)
  error('SLMENGINE:inconsistentdata','x and y must be the same size')
end

% were there any NaN or inf elements in the data?
k = isnan(x) | isnan(y) | isinf(x) | isinf(y);
if any(k)
  % drop them from the analysis
  x(k) = [];
  y(k) = [];
  
  % also drop corresponding weights if they were supplied on inf or NaN
  % elements.
  if ~isempty(prescription.Weights)
    prescription.Weights(k) = [];
  end
  n = length(x);
end

% if weights or errorbars were set, verify the sizes of these
% parameters, compared to the number of data points.
if ~isempty(prescription.Weights)
  prescription.Weights = prescription.Weights(:);
  if n~=length(prescription.Weights)
    error('SLMENGINE:inconsistentweights','Weights vector must be the same length as # of data points')
  end
end
if ~isempty(prescription.ErrorBar)
  EB = prescription.ErrorBar;
  if length(EB) == 1
    prescription.ErrorBar = repmat(EB,n,2);
  elseif (size(EB,1) == 1)
    prescription.ErrorBar = repmat(EB',1,2);
  elseif (size(EB,2) == 1)
    prescription.ErrorBar = repmat(EB,1,2);
  end
end

% we need to scale y to minimize any numerical issues.
% note that the scaling does not change the signs of any derivatives,
% so monotonicity and curvature constraints are not accidentally
% inverted.
% scaleproblem modifies y into yhat, so that yhat now lies in the
% interval [1/phi,phi], where phi is the golden ratio. This
% transformation has the property that all dependent values vary
% over a range of 1.
[yhat,prescriptionhat,YScale,YShift] = scaleproblem(x,y,prescription);

% note that we will report the precription here, as well as the
% adjusted one after the y shifting and y scaling has been applied.
if prescription.Verbosity > 1
  disp('=========================================')
  disp('Model Prescription')
  disp(prescription)
  
  if strcmp(prescription.Scaling,'on')
    disp('=========================================')
    disp('Post-scaling Model Prescription')
    disp(prescriptionhat)
  end
end
% we will be using prescriptionhat from here on, because that has any
% shift/scale parameters built into everything.

% is this a free knot problem?
if strcmp(prescriptionhat.InteriorKnots,'free')
  % allow the interior knots to vary, using fmincon
  
  % get the starting values for the knots
  % knots vector, dx
  if length(prescriptionhat.Knots)==1
    [knots,nk] = chooseknots(prescriptionhat.Knots,n,x);
  else
    % we should check that the knots contain the data
    knots = sort(prescriptionhat.Knots(:));
    nk= length(knots);
    if (knots(1)>min(x)) || (knots(end)<max(x))
      error('SLMENGINE:inadequateknots',['Knots do not contain the data. Data range: ',num2str([min(x),max(x)])])
    end
  end
  % there must be at least 3 knots
  if nk<3
    error('SLMENGINE:inadequateknots','Free knot estimation requires at least 3 knots')
  end
  
  % set the constraints for fmincon:
  %   no pair of knots may lie too close to each other
  tol = min(0.001,0.1/(nk-1));
  mindelta = tol*(knots(end) - knots(1));
  A = zeros(2,nk-2);
  % and we want the knots to be in increasing order, but still lie inside
  % the first and final knot.
  A(1,1) = -1;
  A(2,nk-2) = 1;
  b = [-knots(1) - mindelta ; knots(end) - mindelta];
  A = [A;full(spdiags(repmat([1 -1],nk-3,1),[0 1],nk-3,nk-2))];
  b = [b;repmat(-mindelta,nk-3,1)];
  
  % set up the optimization parameters (fmincon)
  fminconoptions = optimset('fmincon');
  fminconoptions.LargeScale = 'off';
  fminconoptions.Algorithm = prescription.FminconAlgorithm;
  if prescription.Verbosity>1
    fminconoptions.Display = 'iter';
  else
    fminconoptions.Display = 'off';
  end
  
  % change the prescription for the subsequent internal calls
  prescrip = prescriptionhat;
  prescrip.Knots = knots;
  prescrip.InteriorKnots = 'fixed';
  prescrip.Plot = 'off';
  
  % call the optimizer
  intknots = knots(2:(end-1));
  intknots = fmincon(@free_knot_obj,intknots,A,b, ...
    [],[],[],[],[],fminconoptions,x,y,prescrip);
  
  % fit the curve one last time with the final knots
  prescrip.Knots(2:(end-1)) = intknots;
  slm = slmengine(x,yhat,prescrip);
  
else
  % the knots are fixed, just estimate the model

  % do the appropriate fit. break them down into special
  % cases, purely for simplicity of the code.
  switch prescriptionhat.Degree
    case {0 'constant'}
      slm = slmengine_constant(x,yhat,prescriptionhat);
    case {1 'linear'}
      slm = slmengine_linear(x,yhat,prescriptionhat);
    case {3 'cubic'}
slm = slmengine_cubic(x,yhat,prescriptionhat);
  end
  
end

% record the shift and scale parameters as used
slm.stats.YShift = YShift;
slm.stats.YScale = YScale;

% back out any shift & scale factors from the model
if strcmp(prescription.Scaling,'on')
  slm.coef(:,1) = (slm.coef(:,1) - YShift)/YScale;
  if size(slm.coef,2) > 1
    slm.coef(:,2) = slm.coef(:,2)/YScale;
  end
end

% do we need to report anything? I.e., what was the
% prescription.Verbosity setting?
if prescription.Verbosity>0
  disp('=========================================')
  disp('MODEL STATISTICS REPORT')
  disp(['Number of data points:      ',num2str(n)])
  disp(['Scale factor applied to y   ',num2str(YScale)])
  disp(['Shift applied to y          ',num2str(YShift)])
  disp(['Total degrees of freedom:   ',num2str(slm.stats.TotalDoF)])
  disp(['Net degrees of freedom:     ',num2str(slm.stats.NetDoF)])
  disp(['R-squared:                  ',num2str(slm.stats.R2)])
  disp(['Adjusted R-squared:         ',num2str(slm.stats.R2Adj)])
  disp(['RMSE:                       ',num2str(slm.stats.RMSE)])
  disp(['Range of prediction errors: ',num2str(slm.stats.ErrorRange)])
  disp(['Error quartiles (25%, 75%): ',num2str(slm.stats.Quartiles)])
  disp('=========================================')
end

% append the shape prescription to the final structure
% this can be used as a template for fitting other functions,
% but also as documentation for the curve fit. This is why
% it is attached to the model structure, rather than as a
% second return argument.
slm.prescription = prescription;
% also attach the data used to build the model, purely
% for documentation purposes.
slm.x = x;
slm.y = y;

% set the extrapolation flag for use by slmeval
slm.Extrapolation = prescription.Extrapolation;

% do we need to plot the curve? I.e., what was the
% prescription.Plot setting?
if strcmp(prescription.Plot,'on')
  % plot the curve
  plotslm(slm)
end

% do we generate points on the curve?
if isempty(prescription.Predictions)
  % none were asked for 
  xp = [];
  yp = [];
else
  % generate the required points on the curve
  xp = linspace(min(x),max(x),prescription.Predictions);
  yp = slmeval(xp,slm,0);
end

% does the user want the result in a 'pp' form?
if strcmpi(prescription.Result,'pp') && ~strcmpi(slm.form,'pp')
  % convert to pp form
  stats = slm.stats;
  slm = slm2pp(slm);
  
  % make sure the stats field gets copied into the pp form
  slm.stats = stats;
end



% ========================================================
% =========== slmengines for each model degree ============
% ========================================================

% ========================================================
% ============== piecewise constant model ================
% ========================================================
function slm = slmengine_constant(x,y,prescription)
% fits a piecewise constant shape prescriptive model

% check for inappropriate properties for a piecewise
% constant model
property_check(prescription,'constant')

% simple things about data first...
% slmengine has already made it a column vector,
% and ensured compatibility in length between
% x and y
nx = length(x);
% knots vector, dx
if length(prescription.Knots)==1
  [knots,nk] = chooseknots(prescription.Knots,nx,x);
else
  % we should check that the knots contain the data
  knots = sort(prescription.Knots(:));
  nk= length(knots);
  if (knots(1)>min(x)) || (knots(end)<max(x))
    error('SLMENGINE:inadequateknots',['Knots do not contain the data. Data range: ',num2str([min(x),max(x)])])
  end
end
dx = diff(knots);
if any(dx==0)
  error('SLMENGINE:indistinctknots','Knots must be distinct.')
end

% number of coefficients to estimate
% a piecewise constant function has one coefficient
% at each knot, but the last knot is "not".
nc = nk-1;

% create empty design, equality, inequality
% constraints arrays. also rhs vectors for each
Mdes = zeros(0,nc);
Meq = Mdes;
Mineq = Mdes;
Mreg = Mdes;

rhseq = [];
rhsineq = [];
rhsreg = [];

% -------------------------------------
% build design matrix - 
% first, bin the data - histc wll do it
[junk,xbin] = histc(x,knots); %#ok
% any point which falls at the top end, is said to
% be in the last bin.
xbin(xbin==nk)=nk-1;
% the design matrix is easy to build - one call to sparse
Mdes = accumarray([(1:nx)',xbin],1,[nx,nc]);
rhs = y;

% have we been given a constraint on the sum of the residuals?
if ~isempty(prescription.SumResiduals)
  Meq = [Meq;sum(Mdes,1)];
  rhseq = [rhseq;prescription.SumResiduals + sum(rhs)];
end

% apply weights
W = prescription.Weights;
if ~isempty(W)
  W = W(:);
  if length(W)~=nx
    error('SLMENGINE:inconsistentweights','Weight vector is not the same length as data')
  end
  
  % scale weight vector
  W = sqrt(nx)*W./norm(W);
  
  W = spdiags(W,0,nx,nx);
  Mdes = W*Mdes;
  rhs = W*rhs;
end

% -------------------------------------
% end conditions do not apply for a piecewise
% constant function

% -------------------------------------
% build regularizer
if nc>2
  % use a simple second order finite difference
  dx1 = dx(1:(nc-2));
  dx2 = dx(2:(nc-1));
  fda = [-2./(dx1.*(dx1+dx2)), 2./(dx1.*dx2), ...
         -2./(dx2.*(dx1+dx2))];
  
  Mreg = zeros(nc-2,nc);
  for i=1:(nc-2);
    Mreg(i,i+[0 1 2]) = fda(i,:);
  end
  rhsreg = zeros(nc-2,1);
end
% scale the regularizer before we apply the
% regularization parameter.
Mreg = Mreg/norm(Mreg,1);

% -------------------------------------
% single point equality constraints
% left hand side
if ~isempty(prescription.LeftValue)
  M = zeros(1,nc);
  M(1) = 1;
  Meq = [Meq;M];
  rhseq = [rhseq;prescription.LeftValue];
end

% right hand side
if ~isempty(prescription.RightValue)
  M = zeros(1,nc);
  M(end) = 1;
  Meq = [Meq;M];
  rhseq = [rhseq;prescription.RightValue];
end

% force curve through an x-y pair
xy = prescription.XY;
if ~isempty(xy)
  n = size(xy,1);
  if any(xy(:,1)<knots(1)) || any(xy(:,1)>knots(end))
    error('SLMENGINE:improperconstraint','XY pairs to force the curve through must lie inside the knots')
  end
  
  [junk,ind] = histc(xy(:,1),knots); %#ok
  ind(ind==(nc+1))=nc;
  M = sparse((1:n)',ind,1,n,nc);
  Meq = [Meq;M];
  rhseq = [rhseq;xy(:,2)];
end

% -------------------------------------
% Integral equality constraint
if ~isempty(prescription.Integral)
  % Rectangle rule. The last knot point
  % adds no contribution to the area, but
  % there are only nk-1 knots anyway.
  M = dx(:)';
  Meq = [Meq;M];
  rhseq = [rhseq;prescription.Integral];
end

% -------------------------------------
% single point inequality constraints
% left hand side minimum
if ~isempty(prescription.LeftMinValue)
  M = zeros(1,nc);
  M(1) = -1;
  Mineq = [Mineq;M];
  rhsineq = [rhsineq;-prescription.LeftMinValue];
end
% left hand side maximum
if ~isempty(prescription.LeftMaxValue)
  M = zeros(1,nc);
  M(1) = 1;
  Mineq = [Mineq;M];
  rhsineq = [rhsineq;prescription.LeftMinValue];
end

% right hand side min
if ~isempty(prescription.RightMinValue)
  M = zeros(1,nc);
  M(end) = -1;
  Mineq = [Mineq;M];
  rhsineq = [rhsineq;-prescription.RightMinValue];
end

% right hand side max
if ~isempty(prescription.RightMaxValue)
  M = zeros(1,nc);
  M(end) = 1;
  Mineq = [Mineq;M];
  rhsineq = [rhsineq;prescription.RightMaxValue];
end

% -------------------------------------
% Error bar inequality constraints
if ~isempty(prescription.ErrorBar)
  % lower bounds
  Mineq = [Mineq;-Mdes];
  rhsineq = [rhsineq;-(y-prescription.ErrorBar(:,1))];
  
  % upper bounds
  Mineq = [Mineq;Mdes];
  rhsineq = [rhsineq;y+prescription.ErrorBar(:,2)];
end

% -------------------------------------
% constant region(s)
if ~isempty(prescription.ConstantRegion)
  % there was at least one constant regions specified
  cr = prescription.ConstantRegion;
  M = zeros(0,nc);
  n = 0;
  for i=1:size(cr,1)
    % enforce constancy between the given range limits
    for j=1:(nc-1)
      if (knots(j+1)>=cr(i,1)) && (knots(j+1)<cr(i,2))
        n=n+1;
        M(n,j+[0 1]) = [-1 1];
      end
    end
  end
  
  Meq = [Meq;M];
  rhseq = [rhseq;zeros(size(M,1),1)];
  
end

% -------------------------------------
% overall inequalities
% global min value
if ~isempty(prescription.MinValue)
  M = -eye(nc,nc);
  Mineq = [Mineq;M];
  rhsineq = [rhsineq;repmat(-prescription.MinValue,nc,1)];
end

% global max value
if ~isempty(prescription.MaxValue)
  M = eye(nc,nc);
  Mineq = [Mineq;M];
  rhsineq = [rhsineq;repmat(prescription.MaxValue,nc,1)];
end

% SimplePeak and SimpleValley are really just composite
% properties. A peak at x == a is equivalent to a monotone
% increasing function for x<=a, and a monotone decreasing
% function for x>=a. Likewise a valley is just the opposite.
% 
% So specifying a peak or a valley will cause any monotonicity
% specs to be ignored. slmset has already checked for this,
% and turned off any other monotonicity based properties.
sp = prescription.SimplePeak;
if isnumeric(sp) && ~isempty(sp)
  prescription.Increasing = [knots(1),sp];
  prescription.Decreasing = [sp,knots(end)];
end
sv = prescription.SimpleValley;
if isnumeric(sv) && ~isempty(sv)
  prescription.Decreasing = [knots(1), sv];
  prescription.Increasing = [sv,knots(end)];
end

% monotonicity?
% increasing regions
incR = prescription.Increasing;
L=0;
if ischar(incR)
  if strcmp(incR,'on')
    L=L+1;
    mono(L).knotlist = 1:nc;
    mono(L).direction = 1;
    mono(L).range = [];
  end
elseif ~isempty(incR)
  for i=1:size(incR,1)
    L=L+1;
    mono(L).knotlist = []; %#ok
    mono(L).direction = 1; %#ok
    mono(L).range = sort(incR(i,:)); %#ok
  end
end
% decreasing regions
decR = prescription.Decreasing;
if ischar(decR)
  if strcmp(decR,'on')
    L=L+1;
    mono(L).knotlist = 1:nc;
    mono(L).direction = -1;
    mono(L).range = [];
  end
elseif ~isempty(decR)
  for i=1:size(decR,1)
    L=L+1;
    mono(L).knotlist = [];
    mono(L).direction = -1;
    mono(L).range = sort(decR(i,:));
  end
end
if L>0
  % there were at least some monotone regions specified
  M = zeros(0,nc);
  n = 0;
  for i=1:L
    if isempty(mono(L).range)
      % the entire range was specified to be monotone
      for j=1:(nc-1)
        n=n+1;
        M(n,j+[0 1]) = [1 -1]*mono(i).direction;
      end
    else
      % only enforce monotonicity between the given range limits
      for j=1:(nc-1)
        if (knots(j)<mono(i).range(2)) && ...
            (knots(j+1)>=mono(i).range(1))
          n=n+1;
          M(n,j+[0 1]) = [1 -1]*mono(i).direction;
        end
      end
      
    end
  end
  
  Mineq = [Mineq;M];
  rhsineq = [rhsineq;zeros(size(M,1),1)];
  
end

% -------------------------------------
% Use envelope inequalities?
switch prescription.Envelope
case 'off'
  % nothing to do in this case
case 'supremum'
  % yhat - y >= 0
  Mineq = [Mineq;-Mdes];
  rhsineq = [rhsineq;-rhs];
case 'infimum'
  % yhat - y <= 0
  Mineq = [Mineq;Mdes];
  rhsineq = [rhsineq;rhs];
end

% -------------------------------------
% scale equalities for unit absolute row sum
if ~isempty(Meq)
  rs = diag(1./sum(abs(Meq),2));
  Meq = rs*Meq;
  rhseq = rs*rhseq;
end
% scale inequalities for unit absolute row sum
if ~isempty(Mineq)
  rs = diag(1./sum(abs(Mineq),2));
  Mineq = rs*Mineq;
  rhsineq = rs*rhsineq;
end

% -------------------------------------
% now worry about the regularization. There are three
% possible cases.
% 1. We have a given regularization parameter
% 2. We have a given rmse that we wish to match
% 3. We must use cross validation to choose the parameter
RP = prescription.Regularization;
if (isnumeric(RP) && (RP>=0)) || ((ischar(RP)) && (strcmpi(RP,'smoothest')))
  % solve the problem using the given regularization parameter
  coef = solve_slm_system(RP,Mdes,rhs,Mreg,rhsreg, ...
    Meq,rhseq,Mineq,rhsineq,prescription);
  
  finalRP = RP;
  
elseif isnumeric(RP) && (RP<0)
  % we must match abs(RP) as the rmse.
  aim_rmse = abs(RP);
  fminbndoptions = optimset('fminbnd');
  RP = fminbnd(@match_rmse,-6,6,fminbndoptions, ...
    Mdes,rhs,Mreg,rhsreg,Meq,rhseq,Mineq,rhsineq, ...
    aim_rmse,prescription);
  
  % we logged the parameter in the optimization. undo that for
  % the final call
  finalRP = 10^RP;
  
  % do one final call to get the final coefficients
  coef = solve_slm_system(RP,Mdes,rhs,Mreg,rhsreg, ...
    Meq,rhseq,Mineq,rhsineq,prescription);
  
elseif ischar(RP)
  % its a cross validation problem to solve. drop out each point
  % in turn from the model, then use fminbnd to minimize the
  % predicted sum of squares at the missing points.
  fminbndoptions = optimset('fminbnd');
  % fminbndoptions.Display = 'iter';
  RP = fminbnd(@min_cv,-6,6,fminbndoptions, ...
    Mdes,rhs,Mreg,rhsreg,Meq,rhseq,Mineq,rhsineq,prescription);
  
  % we logged the parameter in the optimization. undo that for
  % the final call
  finalRP = 10^RP;
  
  % do one final call to get the final coefficients
  coef = solve_slm_system(RP,Mdes,rhs,Mreg,rhsreg, ...
    Meq,rhseq,Mineq,rhsineq,prescription);
end

% -------------------------------------
% unpack coefficients into the result structure
slm.form = 'slm';
slm.degree = 0;
slm.knots = knots;
slm.coef = coef;

% degrees of freedom available
slmstats.TotalDoF = nk - 1;
slmstats.NetDoF = slmstats.TotalDoF - size(Meq,1);
% this function does all of the stats, stuffing into slmstats
slm.stats = modelstatistics(slmstats,y,coef,prescription.YScale,Mdes);

slm.stats.finalRP = finalRP;



% ========================================================
% =============== piecewise linear model =================
% ========================================================
function slm = slmengine_linear(x,y,prescription)
% fits a piecewise linear shape prescriptive model

% check for inappropriate properties for a piecewise
% linear model
property_check(prescription,'linear')

% simple things about data first...
% slmengine has already made it a column vector,
% and ensured compatibility in length between
% x and y
nx = length(x);

% knots vector, dx
if length(prescription.Knots)==1
  % just given a number of knots
  [knots,nk] = chooseknots(prescription.Knots,nx,x);
else
  % we should check that the knots contain the data
  knots = sort(prescription.Knots(:));
  nk= length(knots);
  if (knots(1)>min(x)) || (knots(end)<max(x))
    error('SLMENGINE:inadequateknots',['Knots do not contain the data. Data range: ',num2str([min(x),max(x)])])
  end
end
dx = diff(knots);
if any(dx==0)
  error('SLMENGINE:indistinctknots','Knots must be distinct.')
end




% number of coefficients to estimate.
% a piecewise linear Hermite has one coefficient
% at each knot.
nc = nk;

% create empty design, equality, inequality
% constraints arrays. also rhs vectors for each
Mdes = zeros(0,nc);
Meq = Mdes;
Mineq = Mdes;
Mreg = Mdes;

rhseq = [];
rhsineq = [];
rhsreg = [];

% -------------------------------------
% build design matrix - 
% first, bin the data - histc wll do it
[junk,xbin] = histc(x,knots); %#ok
% any point which falls at the top end, is said to
% be in the last bin.
xbin(xbin==nk)=nk-1;

% build design matrix
t = (x - knots(xbin))./dx(xbin);
Mdes = accumarray([repmat((1:nx)',2,1), ...
  [xbin;xbin+1]],[1-t;t],[nx,nc]);
rhs = y;

% have we been given a constraint on the sum of the residuals?
if ~isempty(prescription.SumResiduals)
  Meq = [Meq;sum(Mdes,1)];
  rhseq = [rhseq;prescription.SumResiduals + sum(rhs)];
end

% apply weights
W = prescription.Weights;
if ~isempty(W)
  W = W(:);
  if length(W)~=nx
    error('SLMENGINE:inconsistentweights','Weight vector is not the same length as data')
  end
  
  % scale weight vector
  W = sqrt(nx)*W./norm(W);
  
  W = spdiags(W,0,nx,nx);
  Mdes = W*Mdes;
  rhs = W*rhs;
end

% -------------------------------------
% build regularizer, only bother if at least 3 knots
if nc>2
  % use a simple second order finite difference
  dx1 = dx(1:(nc-2));
  dx2 = dx(2:(nc-1));
  fda = [-2./(dx1.*(dx1+dx2)), 2./(dx1.*dx2), ...
         -2./(dx2.*(dx1+dx2))];
  
  Mreg = zeros(nc-2,nc);
  for i=1:(nc-2);
    Mreg(i,i+[0 1 2]) = fda(i,:);
  end
  rhsreg = zeros(nc-2,1);
end
% scale the regularizer before we apply the
% regularization parameter.
Mreg = Mreg/norm(Mreg,1);

% -------------------------------------
% end conditions only apply for a piecewise linear
% function, IF the conditions are periodicity
if strcmp(prescription.EndConditions,'periodic')
  % set the first and last function values equal
  M = zeros(1,nc);
  M([1,end]) = [-1 1];
  Meq = [Meq;M];
  rhseq = [rhseq;0];
end

% -------------------------------------
% single point equality constraints
% left hand side
if ~isempty(prescription.LeftValue)
  M = zeros(1,nc);
  M(1) = 1;
  Meq = [Meq;M];
  rhseq = [rhseq;prescription.LeftValue];
end

% right hand side
if ~isempty(prescription.RightValue)
  M = zeros(1,nc);
  M(end) = 1;
  Meq = [Meq;M];
  rhseq = [rhseq;prescription.RightValue];
end

% left end slope
if ~isempty(prescription.LeftSlope)
  M = zeros(1,nc);
  M(1:2) = [-1 1]/dx(1);
  Meq = [Meq;M];
  rhseq = [rhseq;prescription.LeftSlope];
end

% Right end slope
if ~isempty(prescription.RightSlope)
  M = zeros(1,nc);
  M([nc-1, nc]) = [-1 1]/dx(end);
  Meq = [Meq;M];
  rhseq = [rhseq;prescription.RightSlope];
end

% force curve through an x-y pair
xy = prescription.XY;
if ~isempty(xy)
  n = size(xy,1);
  if any(xy(:,1)<knots(1)) || any(xy(:,1)>knots(end))
    error('SLMENGINE:improperconstraint','XY pairs to force the curve through must lie inside the knots')
  end
  
  [junk,ind] = histc(xy(:,1),knots); %#ok
  ind(ind==(nc+1))=nc;
  t = (xy(:,1) - knots(ind))./dx(ind);
  
  M = sparse(repmat((1:n)',1,2),[ind,ind+1],[1-t,t],n,nc);
  Meq = [Meq;M];
  rhseq = [rhseq;xy(:,2)];
end

% force slope at a point, or set of points
xyp = prescription.XYP;
if ~isempty(xyp)
  n = size(xyp,1);
  if any(xyp(:,1)<knots(1)) || any(xyp(:,1)>knots(end))
    error('SLMENGINE:improperconstraint','X-Y'' pairs to enforce the slope at must lie inside the knots')
  end
  
  [junk,ind] = histc(xyp(:,1),knots); %#ok
  ind(ind==(nc+1))=nc;
  
  M = sparse(repmat((1:n)',1,2),[ind,ind+1],(1./dx(ind))*[-1,1],n,nc);
  Meq = [Meq;M];
  rhseq = [rhseq;xyp(:,2)];
  
end

% -------------------------------------
% constant region(s)
if ~isempty(prescription.ConstantRegion)
  % there was at least one constant region specified
  cr = prescription.ConstantRegion;
  M = zeros(0,nc);
  n = 0;
  for i=1:size(cr,1)
    % enforce constancy between the given range limits
    for j=1:(nc-1)
      if (knots(j+1)>=cr(i,1)) && (knots(j)<cr(i,2))
        n=n+1;
        M(n,j+[0 1]) = [-1 1];
      end
    end
  end
  
  Meq = [Meq;M];
  rhseq = [rhseq;zeros(size(M,1),1)];
  
end

% -------------------------------------
% linear region(s)
% a linear region implies that across any knot
% inside that region, the slopes must be constant
if ~isempty(prescription.LinearRegion) && (nc>2)
  % there was at least one linear region specified
  lr = prescription.LinearRegion;
  M = zeros(0,nc);
  n = 0;
  for i=1:size(lr,1)
    % enforce constancy of slope between the given range limits
    for j=1:(nc-2)
      if (knots(j+1)>=lr(i,1)) && (knots(j+1)<lr(i,2))
        n=n+1;
        M(n,j+[0 1 2]) = [-1 1 0]/dx(j) - [0 -1 1]/dx(j+1);
      end
    end
  end
  
  Meq = [Meq;M];
  rhseq = [rhseq;zeros(size(M,1),1)];
end

% -------------------------------------
% constant knot intervals(s)
if ~isempty(prescription.SegmentConstant)
  % there was at least one constant knot segment specified
  cindex = prescription.SegmentConstant;
  ni = numel(cindex);
  M = zeros(ni,nc);
  for i=1:ni
    % enforce constancy on a given knot interval or intervals
    M(i,cindex(i)+[0 1]) = [-1 1];
  end
  
  Meq = [Meq;M];
  rhseq = [rhseq;zeros(size(M,1),1)];
end

% -------------------------------------
% Integral equality constraint
if ~isempty(prescription.Integral)
  % Trapezoidal rule over the support.
  M = ([dx(:)',0] + [0,dx(:)'])/2;
  Meq = [Meq;M];
  rhseq = [rhseq;prescription.Integral];
end

% -------------------------------------
% single point inequality constraints
% left hand side minimum
if ~isempty(prescription.LeftMinValue)
  M = zeros(1,nc);
  M(1) = -1;
  Mineq = [Mineq;M];
  rhsineq = [rhsineq;-prescription.LeftMinValue];
end
% left hand side maximum
if ~isempty(prescription.LeftMaxValue)
  M = zeros(1,nc);
  M(1) = 1;
  Mineq = [Mineq;M];
  rhsineq = [rhsineq;prescription.LeftMinValue];
end

% right hand side min
if ~isempty(prescription.RightMinValue)
  M = zeros(1,nc);
  M(end) = -1;
  Mineq = [Mineq;M];
  rhsineq = [rhsineq;-prescription.RightMinValue];
end

% right hand side max
if ~isempty(prescription.RightMaxValue)
  M = zeros(1,nc);
  M(end) = 1;
  Mineq = [Mineq;M];
  rhsineq = [rhsineq;prescription.RightMaxValue];
end

% left end slope min
if ~isempty(prescription.LeftMinSlope)
  M = zeros(1,nc);
  M(1:2) = [1 -1]/dx(1);
  Mineq = [Mineq;M];
  rhsineq = [rhsineq;-prescription.LeftMinSlope];
end

% left end slope max
if ~isempty(prescription.LeftMaxSlope)
  M = zeros(1,nc);
  M(1:2) = [-1 1]/dx(1);
  Mineq = [Mineq;M];
  rhsineq = [rhsineq;prescription.LeftMaxSlope];
end

% right end slope min
if ~isempty(prescription.RightMinSlope)
  M = zeros(1,nc);
  M([nc-1,nc]) = [1 -1]/dx(end);
  Mineq = [Mineq;M];
  rhsineq = [rhsineq;-prescription.RightMinSlope];
end

% right end slope max
if ~isempty(prescription.RightMaxSlope)
  M = zeros(1,nc);
  M([nc-1,nc]) = [-1 1]/dx(end);
  Mineq = [Mineq;M];
  rhsineq = [rhsineq;prescription.RightMaxSlope];
end

% -------------------------------------
% overall inequalities
% global min value
if ~isempty(prescription.MinValue)
  M = -eye(nc,nc);
  Mineq = [Mineq;M];
  rhsineq = [rhsineq;repmat(-prescription.MinValue,nc,1)];
end

% global max value
if ~isempty(prescription.MaxValue)
  M = eye(nc,nc);
  Mineq = [Mineq;M];
  rhsineq = [rhsineq;repmat(prescription.MaxValue,nc,1)];
end

% global min slope
if ~isempty(prescription.MinSlope)
  M = full(spdiags((1./dx)*[1 -1],[0 1],nc-1,nc));
  Mineq = [Mineq;M];
  rhsineq = [rhsineq;repmat(-prescription.MinSlope,nc-1,1)];
end

% global max slope
if ~isempty(prescription.MaxSlope)
  M = full(spdiags((1./dx)*[-1 1],[0 1],nc-1,nc));
  Mineq = [Mineq;M];
  rhsineq = [rhsineq;repmat(prescription.MaxSlope,nc-1,1)];
end

% -------------------------------------
% SimplePeak and SimpleValley are really just composite
% properties. A peak at x == a is equivalent to a monotone
% increasing function for x<=a, and a monotone decreasing
% function for x>=a. Likewise a valley is just the opposite.
% 
% So specifying a peak or a valley will cause any monotonicity
% specs to be ignored. slmset has already checked for this,
% and turned off any other monotonicity based properties.
sp = prescription.SimplePeak;
if isnumeric(sp) && ~isempty(sp)
  prescription.Increasing = [knots(1),sp];
  prescription.Decreasing = [sp,knots(end)];
end
sv = prescription.SimpleValley;
if isnumeric(sv) && ~isempty(sv)
  prescription.Decreasing = [knots(1), sv];
  prescription.Increasing = [sv,knots(end)];
end

% -------------------------------------
% monotonicity?
% increasing regions
incR = prescription.Increasing;
L=0;
if ischar(incR)
  if strcmp(incR,'on')
    L=L+1;
    mono(L).knotlist = 'all';
    mono(L).direction = 1;
    mono(L).range = [];
  end
elseif ~isempty(incR)
  for i=1:size(incR,1)
    L=L+1;
    mono(L).knotlist = []; %#ok
    mono(L).direction = 1; %#ok
    mono(L).range = sort(incR(i,:)); %#ok
  end
end
% decreasing regions
decR = prescription.Decreasing;
if ischar(decR)
  if strcmp(decR,'on')
    L=L+1;
    mono(L).knotlist = 'all';
    mono(L).direction = -1;
    mono(L).range = [];
  end
elseif ~isempty(decR)
  for i=1:size(decR,1)
    L=L+1;
    mono(L).knotlist = [];
    mono(L).direction = -1;
    mono(L).range = sort(decR(i,:));
  end
end
if L>0
  % there were at least some monotone regions specified
  M = zeros(0,nc);
  n = 0;
  for i=1:L
    if isempty(mono(L).range)
      % the entire range was specified to be monotone
      for j=1:(nc-1)
        n=n+1;
        M(n,j+[0 1]) = [1 -1]*mono(i).direction;
      end
    else
      % only enforce monotonicity between the given range limits
      for j=1:(nc-1)
        if (knots(j)<mono(i).range(2)) && ...
            (knots(j+1)>=mono(i).range(1))
          n=n+1;
          M(n,j+[0 1]) = [1 -1]*mono(i).direction;
        end
      end
      
    end
  end
  
  Mineq = [Mineq;M];
  rhsineq = [rhsineq;zeros(size(M,1),1)];
end

% -------------------------------------
% PositiveInflection and NegativeInflection are really just
% composite properties. A point of inflection at x == a is
% equivalent to a change in the sign of the curvature at
% x == a.
pinf = prescription.PositiveInflection;
if isnumeric(pinf) && ~isempty(pinf)
  prescription.ConcaveDown = [knots(1),pinf];
  prescription.ConcaveUp = [pinf,knots(end)];
end
ninf = prescription.NegativeInflection;
if isnumeric(ninf) && ~isempty(ninf)
  prescription.ConcaveUp = [knots(1),ninf];
  prescription.ConcaveDown = [ninf,knots(end)];
end

% -------------------------------------
% concave up or down?
% regions of postive curvature
cuR = prescription.ConcaveUp;
L=0;
if ischar(cuR)
  if strcmp(cuR,'on')
    L=L+1;
    curv(L).knotlist = 'all';
    curv(L).direction = 1;
    curv(L).range = [];
  end
elseif ~isempty(cuR)
  for i=1:size(cuR,1)
    L=L+1;
    curv(L).knotlist = []; %#ok
    curv(L).direction = 1; %#ok
    curv(L).range = sort(cuR(i,:)); %#ok
  end
end

% decreasing regions
cdR = prescription.ConcaveDown;
if ischar(cdR)
  if strcmp(cdR,'on')
    L=L+1;
    curv(L).knotlist = 'all';
    curv(L).direction = -1;
    curv(L).range = [];
  end
elseif ~isempty(cdR)
  for i=1:size(cdR,1)
    L=L+1;
    curv(L).knotlist = [];
    curv(L).direction = -1;
    curv(L).range = sort(cdR(i,:));
  end
end
if L>0
  % there were at least some regions with specified curvature
  M = zeros(0,nc);
  n = 0;
  for i=1:L
    if isempty(curv(L).range)
      % the entire range was specified to be
      % curved in some direction
      for j=2:(nc-1)
        n=n+1;
        M(n,j+[-1 0 1]) = ([-1 1 0]/dx(j-1) - [0 -1 1]/dx(j))*curv(i).direction;
      end
    else
      % only enforce curvature between the given range limits
      for j=2:(nc-1)
        if (knots(j)<curv(i).range(2)) && ...
            (knots(j)>=curv(i).range(1))
          n=n+1;
          M(n,j+[-1 0 1]) = ([-1 1 0]/dx(j-1) - [0 -1 1]/dx(j))*curv(i).direction;
        end
      end
      
    end
  end
  
  Mineq = [Mineq;M];
  rhsineq = [rhsineq;zeros(size(M,1),1)];
end

% -------------------------------------
% Error bar inequality constraints
if ~isempty(prescription.ErrorBar)
  % lower bounds
  Mineq = [Mineq;-Mdes];
  rhsineq = [rhsineq;-(y-prescription.ErrorBar(:,1))];
  
  % upper bounds
  Mineq = [Mineq;Mdes];
  rhsineq = [rhsineq;y+prescription.ErrorBar(:,2)];
end

% -------------------------------------
% Use envelope inequalities?
switch prescription.Envelope
case 'off'
  % nothing to do in this case
case 'supremum'
  % yhat - y >= 0
  Mineq = [Mineq;-Mdes];
  rhsineq = [rhsineq;-rhs];
case 'infimum'
  % yhat - y <= 0
  Mineq = [Mineq;Mdes];
  rhsineq = [rhsineq;rhs];
end

% -------------------------------------
% scale equalities for unit absolute row sum
if ~isempty(Meq)
  rs = diag(1./sum(abs(Meq),2));
  Meq = rs*Meq;
  rhseq = rs*rhseq;
end
% scale inequalities for unit absolute row sum
if ~isempty(Mineq)
  rs = diag(1./sum(abs(Mineq),2));
  Mineq = rs*Mineq;
  rhsineq = rs*rhsineq;
end

% -------------------------------------
% now worry about the regularization. There are four
% possible cases.
% 1. We have a given regularization parameter
% 2. We have a given rmse that we wish to match
% 3. We must use cross validation to choose the parameter
RP = prescription.Regularization;
if (isnumeric(RP) && (RP>=0)) || ((ischar(RP)) && (strcmpi(RP,'smoothest')))
  % solve the problem using the given regularization parameter
  coef = solve_slm_system(RP,Mdes,rhs,Mreg,rhsreg, ...
    Meq,rhseq,Mineq,rhsineq,prescription);
  
  finalRP = RP;
  
elseif isnumeric(RP) && (RP<0)
  % we must match abs(RP) as the rmse.
  aim_rmse = abs(RP);
  fminbndoptions = optimset('fminbnd');
  RP = fminbnd(@match_rmse,-6,6,fminbndoptions, ...
    Mdes,rhs,Mreg,rhsreg,Meq,rhseq,Mineq,rhsineq, ...
    aim_rmse,prescription);
  
  % we logged the parameter in the optimization. undo that for
  % the final call
  finalRP = 10^RP;
  
  % do one final call to get the final coefficients
  coef = solve_slm_system(RP,Mdes,rhs,Mreg,rhsreg, ...
    Meq,rhseq,Mineq,rhsineq,prescription);
  
elseif ischar(RP)
  % its a cross validation problem to solve. drop out each point
  % in turn from the model, then use fminbnd to minimize the
  % predicted sum of squares at the missing points.
  fminbndoptions = optimset('fminbnd');
  % fminbndoptions.Display = 'iter';
  RP = fminbnd(@min_cv,-6,6,fminbndoptions, ...
    Mdes,rhs,Mreg,rhsreg,Meq,rhseq,Mineq,rhsineq,prescription);
  
  % we logged the parameter in the optimization. undo that for
  % the final call
  finalRP = 10^RP;
  
  % do one final call to get the final coefficients
  coef = solve_slm_system(RP,Mdes,rhs,Mreg,rhsreg, ...
    Meq,rhseq,Mineq,rhsineq,prescription);
end

% -------------------------------------
% unpack coefficients into the result structure
slm.form = 'slm';
slm.degree = 1;
slm.knots = knots;
slm.coef = coef;

% generate model statistics
slmstats.TotalDoF = nk;
slmstats.NetDoF = slmstats.TotalDoF - size(Meq,1);
% this function does all of the stats, stuffing into slmstats
slm.stats = modelstatistics(slmstats,y,coef,prescription.YScale,Mdes);

slm.stats.finalRP = finalRP;



% ========================================================
% =============== piecewise cubic model ==================
% ========================================================
function slm = slmengine_cubic(x,y,prescription)
% fits a piecewise cubic shape prescriptive model

% check for inappropriate properties for a cubic model
property_check(prescription,'cubic')

% simple things about data first...
% slmengine has already made it a column vector,
% and ensured compatibility in length between
% x and y
nx = length(x);

% knots vector, dx
if length(prescription.Knots)==1
  % just given a number of knots
  [knots,nk] = chooseknots(prescription.Knots,nx,x);
else
  % we should check that the knots contain the data
  knots = sort(prescription.Knots(:));
  nk= length(knots);
  if (knots(1)>min(x)) || (knots(end)<max(x))
    error('SLMENGINE:inadequateknots',['Knots do not contain the data. Data range: ',num2str([min(x),max(x)])])
  end
end
dx = diff(knots);
if any(dx==0)
  error('SLMENGINE:indistinctknots','Knots must be distinct.')
end

% number of coefficients to estimate.
% a piecewise cubic Hermite has two coefficients
% at each knot.
nc = 2*nk;

% create empty inequality constraints arrays. also
% a rhs vector
Mineq = zeros(0,nc);
rhsineq = [];
Meq = zeros(0,nc);
rhseq = [];

% -------------------------------------
% build design matrix - 
% first, bin the data - histc wll do it best
[junk,xbin] = histc(x,knots); %#ok
% any point which falls at the top end, is said to
% be in the last bin.
xbin(xbin==nk)=nk-1;

% build design matrix
t = (x - knots(xbin))./dx(xbin);
t2 = t.^2;
t3 = t.^3;
s2 = (1-t).^2;
s3 = (1-t).^3;

vals = [3*s2-2*s3 ; 3*t2-2*t3 ; ...
  -(s3-s2).*dx(xbin) ; (t3-t2).*dx(xbin)];

% the coefficients will be stored in two blocks,
% first nk function values, then nk derivatives.
Mdes = accumarray([repmat((1:nx)',4,1), ...
  [xbin;xbin+1;nk+xbin;nk+xbin+1]],vals,[nx,nc]);
rhs = y;

% have we been given a constraint on the sum of the residuals?
if ~isempty(prescription.SumResiduals)
  Meq = [Meq;sum(Mdes,1)];
  rhseq = [rhseq;prescription.SumResiduals + sum(rhs)];
end

% apply weights
W = prescription.Weights;
if ~isempty(W)
  W = W(:);
  if length(W)~=nx
    error('SLMENGINE:inconsistentweights','Weight vector is not the same length as data')
  end
  
  % scale weight vector
  W = sqrt(nx)*W./norm(W);
  
  W = spdiags(W,0,nx,nx);
  Mdes = W*Mdes;
  rhs = W*rhs;
end

% -------------------------------------
% Regularizer
% We are integrating the piecewise linear f''(x), as
% a quadratic form in terms of the (unknown) second
% derivatives at the knots.
Mreg=zeros(nk,nk);
Mreg(1,1:2)=[dx(1)/3 , dx(1)/6];
Mreg(nk,nk+[-1 0])=[dx(end)/6 , dx(end)/3];
for i=2:(nk-1)
  Mreg(i,i+[-1 0 1])=[dx(i-1)/6 , (dx(i-1)+dx(i))/3 , dx(i)/6];
end
% do a matrix square root. cholesky is simplest & ok since regmat is
% positive definite. this way we can write the quadratic form as:
%    s'*r'*r*s,
% where s is the vector of second derivatives at the knots.
Mreg=chol(Mreg);

% next, write the second derivatives as a function of the
% function values and first derivatives.   s = [sf,sd]*[f;d]
sf = zeros(nk,nk);
sd = zeros(nk,nk);
for i = 1:(nk-1)
  sf(i,i+[0 1]) = [-1 1].*(6/(dx(i)^2));
  sd(i,i+[0 1]) = [-4 -2]./dx(i);
end
sf(nk,nk+[-1 0]) = [1 -1].*(6/(dx(end)^2));
sd(nk,nk+[-1 0]) = [2 4]./dx(end);
Mreg = Mreg*[sf,sd];

% scale the regularizer before we apply the
% regularization parameter.
Mreg = Mreg/norm(Mreg,1);
rhsreg = zeros(nk,1);

% -------------------------------------
% C2 continuity across knots
if strcmp(prescription.C2,'on')
  MC2 = zeros(nk-2,nc);
  for i = 1:(nk-2)
    MC2(i,[i,i+1]) = [6 -6]./(dx(i).^2);
    MC2(i,[i+1,i+2]) = MC2(i,[i+1,i+2]) + [6 -6]./(dx(i+1).^2);
    
    MC2(i,nk+[i,i+1]) = [2 4]./dx(i);
    MC2(i,nk+[i+1,i+2]) = MC2(i,nk+[i+1,i+2]) + [4 2]./dx(i+1);
  end
  
  Meq = [Meq;MC2];
  rhseq = [rhseq;zeros(nk-2,1)];
end

% -------------------------------------
% end conditions (do nothing here if they are 'estimate')
if strcmp(prescription.EndConditions,'periodic')
  % set the first and last function values equal
  M = zeros(3,nc);
  M(1,[1,nk]) = [-1 1];
  % as well as the first derivatives
  M(2,[nk+1,2*nk]) = [-1 1];
  % and the second derivatives
  M(3,[1 2]) = [-6 6]/(dx(1).^2);
  M(3,[nk+1,nk+2]) = [-4 -2]/dx(1);
  M(3,[nk-1,nk]) = [-6 6]/(dx(end).^2);
  M(3,[2*nk-1,2*nk]) = [-2 -4]/dx(end);
  
  Meq = [Meq;M];
  rhseq = [rhseq;zeros(3,1)];
elseif strcmp(prescription.EndConditions,'notaknot') && nk==2
  % third derivative is continuous across 2nd knot
  % this is the special case for 2 segments, so 1 internal knot.
  M = zeros(1,nc);
  
  M(1,[1 2]) = [12 -12]/(dx(1).^3);
  M(1,[2 3]) = M(1,2:3) - [12 -12]/(dx(2).^3);
  
  M(1,nk+[1 2]) = [6 6]/(dx(1).^2);
  M(1,nk+[2 3]) = M(1,nk+[2 3]) - [6 6]/(dx(2).^2);
  
  Meq = [Meq;M];
  rhseq = [rhseq;0];
elseif strcmp(prescription.EndConditions,'notaknot') && nk>2
  % third derivative is continuous at 2nd & penultimate knots
  M = zeros(2,nc);
  
  M(1,[1 2]) = [12 -12]/(dx(1).^3);
  M(1,[2 3]) = M(1,2:3) - [12 -12]/(dx(2).^3);
  M(1,nk+[1 2]) = [6 6]/(dx(1).^2);
  M(1,nk+[2 3]) = M(1,nk+[2 3]) - [6 6]/(dx(2).^2);
  
  M(2,nk + [-2 -1]) = [12 -12]/(dx(end - 1).^3);
  M(2,nk + [-1 0]) = M(2,nk + [-1 0]) - [12 -12]/(dx(end).^3);
  M(2,2*nk+[-2 -1]) = [6 6]/(dx(end-1).^2);
  M(2,2*nk+[-1 0]) = M(2,2*nk+[-1 0]) - [6 6]/(dx(end).^2);
  
  Meq = [Meq;M];
  rhseq = [rhseq;0;0];
elseif strcmp(prescription.EndConditions,'natural')
  % second derivative at each end is zero
  M = zeros(2,nc);
  
  M(1,[1 2]) = [-6 6]/(dx(1).^2);
  M(1,nk+[1 2]) = [-4 -2]/dx(1);
  
  M(2,nk+[-1 0]) = [6 -6]/(dx(end).^2);
  M(2,2*nk+[-1 0]) = [2 4]/dx(end);
  
  Meq = [Meq;M];
  rhseq = [rhseq;0;0];
end

% -------------------------------------
% single point equality constraints at a knot
% left hand side
if ~isempty(prescription.LeftValue)
  M = zeros(1,nc);
  M(1,1) = 1;
  Meq = [Meq;M];
  rhseq = [rhseq;prescription.LeftValue];
end

% right hand side
if ~isempty(prescription.RightValue)
  M = zeros(1,nc);
  M(1,nk) = 1;
  Meq = [Meq;M];
  rhseq = [rhseq;prescription.RightValue];
end

% left end slope
if ~isempty(prescription.LeftSlope)
  M = zeros(1,nc);
  M(1,nk+1) = 1;
  Meq = [Meq;M];
  rhseq = [rhseq;prescription.LeftSlope];
end

% Right end slope
if ~isempty(prescription.RightSlope)
  M = zeros(1,nc);
  M(1,2*nk) = 1;
  Meq = [Meq;M];
  rhseq = [rhseq;prescription.RightSlope];
end

% -------------------------------------
% force curve through an x-y pair or pairs
xy = prescription.XY;
if ~isempty(xy)
  n = size(xy,1);
  if any(xy(:,1)<knots(1)) || any(xy(:,1)>knots(end))
    error('SLMENGINE:improperconstraint','XY pairs to force the curve through must lie inside the knots')
  end
  
  [junk,ind] = histc(xy(:,1),knots); %#ok
  ind(ind==(nk))=nk-1;
  
  t = (xy(:,1) - knots(ind))./dx(ind);
  t2 = t.^2;
  t3 = t.^3;
  s2 = (1-t).^2;
  s3 = (1-t).^3;
  
  vals = [3*s2-2*s3 ; 3*t2-2*t3 ; ...
    -(s3-s2).*dx(ind) ; (t3-t2).*dx(ind)];
  
  M = accumarray([repmat((1:n)',4,1), ...
    [ind;ind+1;nk+ind;nk+ind+1]],vals,[n,nc]);
  
  Meq = [Meq;M];
  rhseq = [rhseq;xy(:,2)];
end

% -------------------------------------
% force slope at a point, or a set of points
xyp = prescription.XYP;
if ~isempty(xyp)
  n = size(xyp,1);
  if any(xyp(:,1)<knots(1)) || any(xyp(:,1)>knots(end))
    error('SLMENGINE:improperconstraint','X-Y'' pairs to enforce the slope at must lie inside the knots')
  end
  
  [junk,ind] = histc(xyp(:,1),knots); %#ok
  ind(ind==(nk))=nk-1;
  
  t = (xyp(:,1) - knots(ind))./dx(ind);
  t2 = t.^2;
  s2 = (1-t).^2;
  s = (1-t);
  
  vals = [(-6*s+6*s2)./dx(ind) ; ...
    (6*t-6*t2)./dx(ind) ; (3*s2-2*s) ; (3*t2-2*t)];
  
  M = accumarray([repmat((1:n)',4,1), ...
    [ind;ind+1;nk+ind;nk+ind+1]],vals,[n,nc]);
  
  Meq = [Meq;M];
  rhseq = [rhseq;xyp(:,2)];
end

% -------------------------------------
% force second derivative at a point, or set of points
xypp = prescription.XYPP;
if ~isempty(xypp)
  n = size(xypp,1);
  if any(xypp(:,1)<knots(1)) || any(xypp(:,1)>knots(end))
    error('SLMENGINE:improperconstraint','X-Y'' pairs to enforce y" at must lie inside the knots')
  end
  
  [junk,ind] = histc(xypp(:,1),knots); %#ok
  ind(ind==(nk))=nk-1;
  
  t = (xypp(:,1) - knots(ind))./dx(ind);
  s = (1-t);
  
  vals = [(6 - 12*s)./dx(ind).^2; ...
    (6 - 12*t)./dx(ind).^2 ; ...
    -(6*s - 2)./dx(ind) ; ...
    (6*t - 2)./dx(ind)];
  
  M = accumarray([repmat((1:n)',4,1), ...
    [ind;ind+1;nk+ind;nk+ind+1]],vals,[n,nc]);
  
  Meq = [Meq;M];
  rhseq = [rhseq;xypp(:,2)];
end

% -------------------------------------
% force third derivative at a point, or set of points
xyppp = prescription.XYPPP;
if ~isempty(xyppp)
  n = size(xyppp,1);
  if any(xyppp(:,1)<knots(1)) || any(xyppp(:,1)>knots(end))
    error('SLMENGINE:improperconstraint','X-Y'''''' pairs to enforce y'''''' at must lie inside the knots')
  end
  
  [junk,ind] = histc(xyppp(:,1),knots); %#ok
  ind(ind==(nk))=nk-1;
  
  vals = [12./dx(ind).^3; -12./dx(ind).^3 ; ...
    6./dx(ind).^2 ;6./dx(ind).^2];
  
  M = accumarray([repmat((1:n)',4,1), ...
    [ind;ind+1;nk+ind;nk+ind+1]],vals,[n,nc]);
  
  Meq = [Meq;M];
  rhseq = [rhseq;xyppp(:,2)];
end

% -------------------------------------
% constant region(s)
if ~isempty(prescription.ConstantRegion)
  % there was at least one constant region specified
  cr = prescription.ConstantRegion;
  M = zeros(0,nc);
  n = 0;
  for i=1:size(cr,1)
    % enforce constancy between the given range limits
    flag = true;
    for j=1:(nk-1)
      if (knots(j+1)>=cr(i,1)) && (knots(j)<cr(i,2))
        % f(j) == f(j+1)
        n=n+1;
        M(n,j+[0 1]) = [-1 1];
        
        % also zero f'' on this interval
        if flag
          % only set this constraint for the first knot
          % interval in a constant region
          n=n+1;
          M(n,j+[0 1 nk nk+1]) = [-6/(dx(j).^2) ; ...
            6/(dx(j).^2) ; -4/dx(j) ; -2/dx(j)];
          
          flag = false;
        end
        
        n=n+1;
        M(n,j+[0 1 nk nk+1]) = [-6/(dx(j).^2) ; ...
          6/(dx(j).^2) ; 2/dx(j) ; 4/dx(j)];
        
      end
    end
  end
  
  Meq = [Meq;M];
  rhseq = [rhseq;zeros(size(M,1),1)];
  
end

% -------------------------------------
% linear region(s)
% a linear region simply means that across any knot
% inside that region, the slopes must be constant
if ~isempty(prescription.LinearRegion) && (nc>2)
  % there was at least one linear region specified
  lr = prescription.LinearRegion;
  M = zeros(0,nc);
  n = 0;
  for i=1:size(lr,1)
    % enforce constancy of slope between the given range limits
    flag = true;
    for j=1:(nk-1)
      if (knots(j+1)>=lr(i,1)) && (knots(j)<lr(i,2))
        % f'(j) == f'(j+1)
        n=n+1;
        M(n,j+nk+[0 1]) = [-1 1];
        
        % also zero f'' at each end of this interval
        if flag
          % only set this constraint for the first knot
          % interval in a linear region
          n=n+1;
          M(n,j+[0 1 nk nk+1]) = [-6/(dx(j).^2) ; ...
            6/(dx(j).^2) ; -4/dx(j) ; -2/dx(j)];
          
          flag = false;
        end
        n=n+1;
        M(n,j+[0 1 nk nk+1]) = [6/(dx(j).^2) ; ...
          -6/(dx(j).^2) ; 2/dx(j) ; 4/dx(j)];
        
      end
    end
  end
  
  Meq = [Meq;M];
  rhseq = [rhseq;zeros(size(M,1),1)];
end

% -------------------------------------
% constant knot intervals(s)
if ~isempty(prescription.SegmentConstant)
  % there was at least one constant knot segment specified
  cindex = prescription.SegmentConstant;
  ni = numel(cindex);
  M = zeros(3*ni,nc);
  for i=1:ni
    % enforce constancy on a given knot interval or intervals
    M(3*i-2,cindex(i) + [0 1]) = [-1 1];
    M(3*i-1,cindex(i) + nk) = 1;
    M(3*i  ,cindex(i) + 1 + nk) = 1;
  end
  
  Meq = [Meq;M];
  rhseq = [rhseq;zeros(size(M,1),1)];
end

% -------------------------------------
% Linear knot intervals(s)
if ~isempty(prescription.SegmentLinear)
  % there was at least one linear knot segment specified
  cindex = prescription.SegmentLinear;
  ni = numel(cindex);
  M = zeros(2*ni,nc);
  for i=1:ni
    % enforce linearity on a given knot interval or intervals
    M(2*i-1,cindex(i) + [0 1]) = [-1 1];
    M(2*i,cindex(i) + nk) = 1;
    
    ci = cindex(i);
    dxi = dx(ci);
    M(2*i-1,ci + [0 1 nk nk+1]) = [-6/(dxi.^2) ; ...
      6/(dxi.^2) ; -4/dxi ; -2/dxi];
    
    M(2*i,ci + [0 1 nk nk+1]) = [6/(dxi.^2) ; ...
      -6/(dxi.^2) ; 2/dxi ; 4/dxi];
  end
  
  Meq = [Meq;M];
  rhseq = [rhseq;zeros(size(M,1),1)];
end

% -------------------------------------
% Integral equality constraint
if ~isempty(prescription.Integral)
  % Simpson's rule over the support will be exact,
  % evaluating at each midpoint. Don't forget that
  % the effective "stepsize" is actually dx(i)/2
  M = zeros(1,nc);
  for i = 1:(nk-1)
    M(1,i+[0 1]) = M(1,i+[0 1]) + dx(i)/6;
    % the midpoint
    M(1,i+[0 1 nk nk+1]) = M(1,i+[0 1 nk nk+1]) + ...
      [1/2 , 1/2 , (1/8).*dx(i) , (-1/8).*dx(i)]*4*dx(i)/6;
  end
  
  Meq = [Meq;M];
  rhseq = [rhseq;prescription.Integral];
end

% -------------------------------------
% single point inequality constraints
% left hand side minimum
if ~isempty(prescription.LeftMinValue)
  M = zeros(1,nc);
  M(1,1) = -1;
  Mineq = [Mineq;M];
  rhsineq = [rhsineq;-prescription.LeftMinValue];
end

% left hand side maximum
if ~isempty(prescription.LeftMaxValue)
  M = zeros(1,nc);
  M(1,1) = 1;
  Mineq = [Mineq;M];
  rhsineq = [rhsineq;prescription.LeftMaxValue];
end

% right hand side min
if ~isempty(prescription.RightMinValue)
  M = zeros(1,nc);
  M(1,nk) = -1;
  Mineq = [Mineq;M];
  rhsineq = [rhsineq;-prescription.RightMinValue];
end

% right hand side max
if ~isempty(prescription.RightMaxValue)
  M = zeros(1,nc);
  M(1,nk) = 1;
  Mineq = [Mineq;M];
  rhsineq = [rhsineq;prescription.RightMaxValue];
end

% left end slope min
if ~isempty(prescription.LeftMinSlope)
  M = zeros(1,nc);
  M(1,nk+1) = -1;
  Mineq = [Mineq;M];
  rhsineq = [rhsineq;-prescription.LeftMinSlope];
end

% left end slope max
if ~isempty(prescription.LeftMaxSlope)
  M = zeros(1,nc);
  M(1,nk+1) = 1;
  Mineq = [Mineq;M];
  rhsineq = [rhsineq;prescription.LeftMaxSlope];
end

% right end slope min
if ~isempty(prescription.RightMinSlope)
  M = zeros(1,nc);
  M(1,nc) = -1;
  Mineq = [Mineq;M];
  rhsineq = [rhsineq;-prescription.RightMinSlope];
end

% right end slope max
if ~isempty(prescription.RightMaxSlope)
  M = zeros(1,nc);
  M(1,nc) = 1;
  Mineq = [Mineq;M];
  rhsineq = [rhsineq;prescription.RightMaxSlope];
end

% -------------------------------------
% overall inequalities
% global min & max value
if ~isempty(prescription.MinValue) || ~isempty(prescription.MaxValue)
  % sample points for necessary conditions are chebychev
  % nodes. No good reason here, except that they are
  % as good a set of points as any other.
  tm = [0.017037; 0.066987; 0.14645; 0.25; 0.37059; ...
    0.5; 0.62941; 0.75; 0.85355; 0.93301; 0.98296];
  nsamp = length(tm);
  
  ntot = nk+(nk-1)*nsamp;
  Mmin = zeros(ntot,nc);
  Mmax = Mmin;
  
  % constrain values at knots
  if ~isempty(prescription.MinValue)
    Mmin(1:nk,1:nk) = -eye(nk,nk);
  end
  if ~isempty(prescription.MaxValue)
    Mmax(1:nk,1:nk) = eye(nk,nk);
  end
  
  % and intermediate sample points
  t2 = tm.^2;
  t3 = tm.^3;
  s2 = (1-tm).^2;
  s3 = (1-tm).^3;
  vals = [3*s2-2*s3 , 3*t2-2*t3 , -s3+s2 , t3-t2];
  
  for j = 1:(nk-1)
    if ~isempty(prescription.MinValue)
      Mmin((1:nsamp) + (j-1)*nsamp + nk , j+[0 1 nk nk+1]) = ...
        -vals*diag([1 1 dx(j) dx(j)]);
    end
    if ~isempty(prescription.MaxValue)
      Mmax((1:nsamp) + (j-1)*nsamp + nk , j+[0 1 nk nk+1]) = ...
        vals*diag([1 1 dx(j) dx(j)]);
    end
  end
  
  if ~isempty(prescription.MinValue)
    Mineq = [Mineq;Mmin];
    rhsineq = [rhsineq;repmat(-prescription.MinValue,ntot,1)];
  end
  if ~isempty(prescription.MaxValue)
    Mineq = [Mineq;Mmax];
    rhsineq = [rhsineq;repmat(prescription.MaxValue,ntot,1)];
  end
end

% global min or max slope
if ~isempty(prescription.MinSlope) || ~isempty(prescription.MaxSlope)
  % its a transformed monotonicity problem. See the
  % explanation provided below for sufficient conditions
  % for monotonicity. The trick here is to build them
  % for a minimum slope. Do that by effectively building
  % an offset into the slope. I.e., assume that theta is the
  % minimum slope allowed. Implicitly subtract the line
  % theta*x from the curve, then use the standard linear
  % system derived below for monotonicity. This will yield
  % a system of linear constraints that are sufficient for
  % a minimum (or maximum) slope of theta.
  
  % 
  MS = zeros(5*(nk - 1) + nk,nc);
  rhsMS = zeros(5*(nk-1) + nk,1);
  
  % constrain the derivatives at each knot
  for n = 1:nk
    MS(n,nk + n) = -1;
  end
  rhsMS(1:nk) = -1;
  
  n = nk + 1;
  for j = 1:(nk - 1)
    % The constraints applied for an increasing function
    % are (in a form that lsqlin will like), i.e., A*x <= b
    %
    %    -delta + theta  <= 0
    %
    % (etc.)
    MS(n,j + [0 1]) = [1 -1];
    rhsMS(n) = -dx(j);
    
    MS(n + 1,j + [0, 1, nk,nk + 1]) = [3, -3, [-1, 1]*dx(j)];
    rhsMS(n + 1) = -3*dx(j);
    
    MS(n + 2,j + [0, 1, nk,nk + 1]) = [3, -3, [1, -1]*dx(j)];
    rhsMS(n + 2) = -3*dx(j);
    
    MS(n + 3,j + [0, 1, nk,nk + 1]) = [9, -9, [1, 2]*dx(j)];
    rhsMS(n + 3) = -6*dx(j);
    
    MS(n + 4,j + [0, 1, nk,nk + 1]) = [9, -9, [2, 1]*dx(j)];
    rhsMS(n + 4) = -6*dx(j);
    
    n = n + 5;
    
  end
  
  % fold these inequalities into any others
  if ~isempty(prescription.MinSlope)
    theta = prescription.MinSlope;
      
    Mineq = [Mineq;MS];
    rhsineq = [rhsineq;theta*rhsMS];
  end
  
  if ~isempty(prescription.MaxSlope)
    % the max slope case is as if we have negated
    % the min slope case.
    theta = prescription.MaxSlope;
    
    Mineq = [Mineq;-MS];
    rhsineq = [rhsineq;-theta*rhsMS];
  end
  
end

% global min & max second derivative
if ~isempty(prescription.MinFpp) || ~isempty(prescription.MaxFpp)
  % just constrain the second derivative over each interval.
  % since it is piecewise linear, that will be both necessary
  % and sufficient. If the curve is C2, then there will be
  % nk inequality constraints needed. If it is only C1, then
  % 2*nk-2.
  isC2 = true;
  if strcmpi(prescription.C2,'off')
    isC2 = false;
  end
  
  if isC2
    % C2, so just constrain the second derivative at each knot point.
    ntot = nk;
    Mmin = zeros(ntot,nc);
    Mmax = Mmin;
    if ~isempty(prescription.MinFpp) 
      Mmin(1,[1, 2, nk+1, nk+2]) = -[[-6 6]/dx(1)^2,[-4 -2]/dx(1)];
      for j = 2:nk
        Mmin(j,j + [-1, 0, nk-1, nk]) = -[[6 -6]/dx(j-1)^2,[2 4]/dx(j-1)];
      end
    end
    if ~isempty(prescription.MaxFpp) 
      Mmax(1,[1, 2, nk+1, nk+2]) = [[-6 6]/dx(1)^2,[-4 -2]/dx(1)];
      for j = 2:nk
        Mmax(j,j + [-1, 0, nk-1, nk]) = [[6 -6]/dx(j-1)^2,[2 4]/dx(j-1)];
      end
    end
    
  else
    % C1, so we need to constrain the second derivative on both sides of
    % each knot in each interval.
    ntot = 2*nk - 2;
    Mmin = zeros(ntot,nc);
    Mmax = Mmin;
    if ~isempty(prescription.MinFpp) 
      for j = 1:(nk-1)
        Mmin(2*j-1,j + [0, 1, nk, nk+1]) = -[[-6 6]/dx(j)^2,[-4 -2]/dx(j)];
        Mmin(2*j,  j + [0, 1, nk, nk+1]) = -[[6 -6]/dx(j)^2,[2 4]/dx(j)];
      end
    end
    if ~isempty(prescription.MaxFpp) 
      for j = 1:(nk-1)
        Mmax(2*j-1,j + [0, 1, nk, nk+1]) = [[-6 6]/dx(j)^2,[-4 -2]/dx(j)];
        Mmax(2*j,  j + [0, 1, nk, nk+1]) = [[6 -6]/dx(j)^2,[2 4]/dx(j)];
      end
    end
    
  end
  
  if ~isempty(prescription.MinFpp)
    Mineq = [Mineq;Mmin];
    rhsineq = [rhsineq;repmat(-prescription.MinFpp,ntot,1)];
  end
  if ~isempty(prescription.MaxFpp)
    Mineq = [Mineq;Mmax];
    rhsineq = [rhsineq;repmat(prescription.MaxFpp,ntot,1)];
  end
end

% -------------------------------------
% SimplePeak and SimpleValley are really just composite
% properties. A peak at x == a is equivalent to a monotone
% increasing function for x<=a, and a monotone decreasing
% function for x>=a. Likewise a valley is just the opposite.
%
% It is best if a peak or valley occurs at a knot.
% 
% So specifying a peak or a valley will cause any monotonicity
% specs to be ignored. slmset has already checked for this,
% and turned off any other monotonicity based properties.
sp = prescription.SimplePeak;
if isnumeric(sp) && ~isempty(sp)
  prescription.Increasing = [knots(1),sp];
  prescription.Decreasing = [sp,knots(end)];
end
sv = prescription.SimpleValley;
if isnumeric(sv) && ~isempty(sv)
  prescription.Decreasing = [knots(1), sv];
  prescription.Increasing = [sv,knots(end)];
end

% -------------------------------------
% monotonicity?
% increasing regions
totalmonotoneintervals = 0;
incR = prescription.Increasing;
L=0;
if ischar(incR)
  if strcmp(incR,'on')
    L=L+1;
    mono(L).knotlist = [1,nk];
    mono(L).direction = 1;
    totalmonotoneintervals = totalmonotoneintervals + nk - 1;
  end
elseif ~isempty(incR)
  for i=1:size(incR,1)
    L=L+1;
    mono(L).knotlist = [max(1,sum(knots <= incR(i,1))), ...
      min(nk,1 + sum(knots < incR(i,2)))]; %#ok
    mono(L).direction = 1; %#ok
    totalmonotoneintervals = totalmonotoneintervals + diff(mono(L).knotlist);
  end
end

% decreasing regions
decR = prescription.Decreasing;
if ischar(decR)
  if strcmp(decR,'on')
    L=L+1;
    mono(L).knotlist = [1,nk];
    mono(L).direction = -1;
    totalmonotoneintervals = totalmonotoneintervals + nk - 1;
  end
elseif ~isempty(decR)
  for i=1:size(decR,1)
    L=L+1;
    mono(L).knotlist = [max(1,sum(knots <= decR(i,1))), ...
      min(nk,1 + sum(knots < decR(i,2)))];
    mono(L).direction = -1;
    totalmonotoneintervals = totalmonotoneintervals + diff(mono(L).knotlist);
  end
end

if L>0
  % there were at least some monotone regions specified
  M = zeros(7*totalmonotoneintervals,nc);
  n = 0;
  for i=1:L
    for j = mono(i).knotlist(1):(mono(i).knotlist(2) - 1)
      % the function must be monotone between
      % knots j and j + 1. The direction over
      % that interval is specified. The constraint
      % system used comes from Fritsch & Carlson, see here:
      %
      % http://en.wikipedia.org/wiki/Monotone_cubic_interpolation
      %
      % Define delta = (y(i+1) - y(i))/(x(i+1) - x(i))
      % Thus delta is the secant slope of the curve across
      % a knot interval. Further, define alpha and beta as
      % the ratio of the derivative at each end of an
      % interval to the secant slope.
      % 
      %  alpha = d(i)/delta
      %  beta = d(i+1)/delta
      %
      % Then we have an elliptically bounded region in the
      % first quadrant that defines the set of monotone cubic
      % segments. We cannot define that elliptical region
      % using a set of linear constraints. However, by use
      % of a system of 7 linear constraints, we can form a
      % set of sufficient conditions such that the curve
      % will be monotone. There will be some few cubic
      % segments that are actually monotone, yet lie outside
      % of the linear system formed. This is acceptable,
      % as our linear approximation here is a sufficient
      % one for monotonicity, although not a necessary one. 
      % It merely says that the spline may be slightly over
      % constrained, i.e., slightly less flexible than is
      % absolutely necessary. (So?)
      %
      % The 7 constraints applied for an increasing function
      % are (in a form that lsqlin will like):
      %
      %    -delta          <= 0
      %    -alpha          <= 0 
      %    -beta           <= 0
      %    -alpha + beta   <= 3
      %     alpha - beta   <= 3
      %     alpha + 2*beta <= 9
      %   2*alpha + beta   <= 9
      %
      % Multiply these inequalities by (y(i+1) - y(i)) to
      % put them into a linear form.
      M(n + 1,j+[0 1]) = [1 -1]*mono(i).direction;
      M(n + 2,nk + j) = -mono(i).direction;
      M(n + 3,nk + j + 1) = -mono(i).direction;
      
      M(n + 4,j + [0, 1, nk,nk + 1]) = mono(i).direction*[3, -3, [-1, 1]*dx(j)];
      M(n + 5,j + [0, 1, nk,nk + 1]) = mono(i).direction*[3, -3, [1, -1]*dx(j)];
      M(n + 6,j + [0, 1, nk,nk + 1]) = mono(i).direction*[9, -9, [1, 2]*dx(j)];
      M(n + 7,j + [0, 1, nk,nk + 1]) = mono(i).direction*[9, -9, [2, 1]*dx(j)];
      
      n = n + 7;
    end
  end
  
  Mineq = [Mineq;M];
  rhsineq = [rhsineq;zeros(size(M,1),1)];
end

% -------------------------------------
% PositiveInflection and NegativeInflection are really just
% composite properties. A point of inflection at x == a is
% equivalent to a change in the sign of the curvature at
% x == a.
pinf = prescription.PositiveInflection;
if isnumeric(pinf) && ~isempty(pinf)
  prescription.ConcaveDown = [knots(1),pinf];
  prescription.ConcaveUp = [pinf,knots(end)];
end
ninf = prescription.NegativeInflection;
if isnumeric(ninf) && ~isempty(ninf)
  prescription.ConcaveUp = [knots(1),ninf];
  prescription.ConcaveDown = [ninf,knots(end)];
end

% -------------------------------------
% concave up or down?
% regions of positive curvature
cuR = prescription.ConcaveUp;
L=0;
if ischar(cuR)
  if strcmp(cuR,'on')
    L=L+1;
    curv(L).knotlist = 'all';
    curv(L).direction = 1;
    curv(L).range = [];
  end
elseif ~isempty(cuR)
  for i=1:size(cuR,1)
    L=L+1;
    curv(L).knotlist = []; %#ok
    curv(L).direction = 1; %#ok
    curv(L).range = sort(cuR(i,:)); %#ok
  end
end
% negative curvature regions
cdR = prescription.ConcaveDown;
if ischar(cdR)
  if strcmp(cdR,'on')
    L=L+1;
    curv(L).knotlist = 'all';
    curv(L).direction = -1;
    curv(L).range = [];
  end
elseif ~isempty(cdR)
  for i=1:size(cdR,1)
    L=L+1;
    curv(L).knotlist = [];
    curv(L).direction = -1;
    curv(L).range = sort(cdR(i,:));
  end
end
if L>0
  % there were at least some regions with specified curvature
  M = zeros(0,nc);
  n = 0;
  for i=1:L
    if isempty(curv(L).range)
      % the entire domain was specified to be
      % curved in some direction
      for j=1:(nk-1)
        n=n+1;
        M(n,j+[0 1]) = curv(i).direction*[6 -6]/(dx(j).^2);
        M(n,nk+j+[0 1]) = curv(i).direction*[4 2]/dx(j);
      end
      n=n+1;
      M(n,nk+[-1 0]) = curv(i).direction*[-6 6]/(dx(end).^2);
      M(n,2*nk+[-1 0]) = curv(i).direction*[-2 -4]/dx(end);
    else
      % only enforce curvature between the given range limits
      % do each knot first.
      for j=1:(nk-1)
        if (knots(j)<curv(i).range(2)) && ...
            (knots(j)>=curv(i).range(1))
          
          n=n+1;
          M(n,j+[0 1]) = curv(i).direction*[6 -6]/(dx(j).^2);
          M(n,nk+j+[0 1]) = curv(i).direction*[4 2]/dx(j);
        end
      end
      
      % also constrain at the endpoints of the range
      curv(i).range = max(min(curv(i).range(:),knots(end)),knots(1));
      [junk,ind] = histc(curv(i).range,knots); %#ok
      ind(ind==(nk))=nk-1;
  
      t = (curv(i).range - knots(ind))./dx(ind);
      s = 1-t;
      
      for j = 1:numel(ind)
        M(n+j,ind(j)+[0 1 nk nk+1]) = -curv(i).direction* ...
          [(6 - 12*s(j))./(dx(ind(j)).^2), (6 - 12*t(j))./(dx(ind(j)).^2) , ...
          (2 - 6*s(j))./dx(ind(j)), (6*t(j) - 2)./dx(ind(j))];
      end
      
      n = n + numel(ind);
    end
  end
  
  Mineq = [Mineq;M];
  rhsineq = [rhsineq;zeros(size(M,1),1)];
end

% -------------------------------------
% Jerk function increasing or decreasing?
% Just a constraint on the sign of the third derivative.
switch prescription.Jerk
  case 'positive'
    fpppsign = 1;
  case 'negative'
    fpppsign = -1;
  otherwise
    fpppsign = 0;
end
if fpppsign
  % constrain the third derivative sign
  M = zeros(nk-1,nc);
  for i=1:(nk-1)
    % one constraint on every interval
    M(i,i+[0 1]) = fpppsign*[-12 12]/(dx(i).^3);
    M(i,nk+i+[0 1]) = fpppsign*[-6 -6]/dx(i).^2;
  end
  
  Mineq = [Mineq;M];
  rhsineq = [rhsineq;zeros(size(M,1),1)];
end

% -------------------------------------
% Error bar inequality constraints
if ~isempty(prescription.ErrorBar)
  % lower bounds
  Mineq = [Mineq;-Mdes];
  rhsineq = [rhsineq;-(y-prescription.ErrorBar(:,1))];
  
  % upper bounds
  Mineq = [Mineq;Mdes];
  rhsineq = [rhsineq;y+prescription.ErrorBar(:,2)];
end

% -------------------------------------
% Use envelope inequalities?
switch prescription.Envelope
case 'off'
  % nothing to do in this case
case 'supremum'
  % yhat - y >= 0
  Mineq = [Mineq;-Mdes];
  rhsineq = [rhsineq;-rhs];
case 'infimum'
  % yhat - y <= 0
  Mineq = [Mineq;Mdes];
  rhsineq = [rhsineq;rhs];
end

% -------------------------------------
% scale equalities for unit absolute row sum
if ~isempty(Meq)
  rs = diag(1./sum(abs(Meq),2));
  Meq = rs*Meq;
  rhseq = rs*rhseq;
end
% scale inequalities for unit absolute row sum
if ~isempty(Mineq)
  rs = diag(1./sum(abs(Mineq),2));
  Mineq = rs*Mineq;
  rhsineq = rs*rhsineq;
end

% -------------------------------------
% now worry about the regularization. There are three
% possible cases.
% 1. We have a given regularization parameter
% 2. We have a given rmse that we wish to match
% 3. We must use cross validation to choose the parameter
RP = prescription.Regularization;
if (isnumeric(RP) && (RP>=0)) || ((ischar(RP)) && (strcmpi(RP,'smoothest')))
  % solve the problem using the given regularization parameter
  
  finalRP = RP;
  [coef,lambda] = solve_slm_system(finalRP,Mdes,rhs,Mreg,rhsreg, ...
    Meq,rhseq,Mineq,rhsineq,prescription);
  
elseif isnumeric(RP) && (RP<0)
  % we must match abs(RP) as the rmse.
  aim_rmse = abs(RP);
  fminbndoptions = optimset('fminbnd');
  RP = fminbnd(@match_rmse,-6,6,fminbndoptions, ...
    Mdes,rhs,Mreg,rhsreg,Meq,rhseq,Mineq,rhsineq, ...
    aim_rmse,prescription);
  
  % we logged the parameter in the optimization. undo that for
  % the final call
  finalRP = 10^RP;
  
  % do one final call to get the final coefficients
  [coef,lambda] = solve_slm_system(finalRP,Mdes,rhs,Mreg,rhsreg, ...
    Meq,rhseq,Mineq,rhsineq,prescription);
  
elseif ischar(RP)
  % its a cross validation problem to solve. drop out each point
  % in turn from the model, then use fminbnd to minimize the
  % predicted sum of squares at the missing points.
  fminbndoptions = optimset('fminbnd');
  % fminbndoptions.Display = 'iter';
  RP = fminbnd(@min_cv,-6,6,fminbndoptions, ...
    Mdes,rhs,Mreg,rhsreg,Meq,rhseq,Mineq,rhsineq,prescription);
  
  % we logged the parameter in the optimization. undo that for
  % the final call
  finalRP = 10^RP;
  
  % do one final call to get the final coefficients
  [coef,lambda] = solve_slm_system(finalRP,Mdes,rhs,Mreg,rhsreg, ...
    Meq,rhseq,Mineq,rhsineq,prescription);
end

% -------------------------------------
% unpack coefficients into the result structure
slm.form = 'slm';
slm.degree = 3;
slm.knots = knots;
slm.coef = reshape(coef,nk,2);

% generate model statistics
slmstats.TotalDoF = 2*nk;
slmstats.NetDoF = slmstats.TotalDoF - size(Meq,1);
% this function does all of the stats, stuffing into slmstats
slm.stats = modelstatistics(slmstats,y,coef,prescription.YScale,Mdes);

slm.stats.finalRP = finalRP;

% ========================================================
% =============== robustfit solution =====================
% ========================================================
function [coef,lambda] = robustfit_slm_system(RP,Mdes,rhs,Mreg,rhsreg,Meq,rhseq,Mineq,rhsineq,prescription)
% solves the final linear system of equations for the model coefficients
% using an iteratively reweighted least squares solution

% we can only get in here if Robust is 'on'. Set it to 'off' now, since
% we don't want recursive calls into this code block.
prescription.Robust = 'off';

% what are the current set of weights? If none, then set to ones.
currentweights = prescription.Weights;

% get a base fit that ignores robustness
[coef,lambda] = solve_slm_system(RP,Mdes,rhs,Mreg,rhsreg,Meq,rhseq,Mineq,rhsineq,prescription);

% compute the residuals
resid = Mdes*coef - rhs;




%{
D = sqrt(eps('double'));
iter = 0;
iterlim = 50;
wxrank = xrank;    % rank of weighted version of x
while((iter==0) || any(abs(b-b0) > D*max(abs(b),abs(b0))))
   iter = iter+1;
   if (iter>iterlim)
      warning(message('stats:statrobustfit:IterationLimit'));
      break;
   end
   
   % Compute residuals from previous fit, then compute scale estimate
   r = y - X*b;
   radj = r .* adjfactor ./ sw;
   s = madsigma(radj,wxrank);
   
   % Compute new weights from these residuals, then re-fit
   w = feval(wfun, radj/(max(s,tiny_s)*tune));
   b0 = b;
   [b(perm),wxrank] = wfit(y,X(:,perm),w);
end
%}














% ========================================================
% ===== (Essentially Stats toolbox code) madsigma ========
% ========================================================
function s = madsigma(r,p)
%MADSIGMA    Compute sigma estimate using MAD of residuals from 0
% used for robust solver
rs = sort(abs(r));
s = median(rs(max(1,p):end)) / 0.6745;

% ========================================================
% ========= aggregate linear system & solve ==============
% ========================================================
function [coef,lambda] = solve_slm_system(RP,Mdes,rhs,Mreg,rhsreg,Meq,rhseq,Mineq,rhsineq,prescription)
% solves the final linear system of equations for the model coefficients 

% trap for the robust solver
if strcmpi(prescription.Robust,'on')
  if prescription.Verbosity > 1
    disp('Robust fit employed')
  end
  [coef,lambda] = robustfit_slm_system(RP,Mdes,rhs,Mreg,rhsreg,Meq,rhseq,Mineq,rhsineq,prescription);
  return
end

if prescription.Verbosity > 1
  % Linear solve parameters
  disp('=========================================')
  disp('LINEAR SYSTEM SOLVER')
  disp(['Design matrix shape:    ',num2str(size(Mdes))])
  disp(['Regularizer shape:      ',num2str(size(Mreg))])
  disp(['Equality constraints:   ',num2str(size(Meq))])
  disp(['Inequality constraints: ',num2str(size(Mineq))])
  disp(['Solver employed:        ',prescription.Solver])
end

% combine design matrix and regularizer
if strcmpi(prescription.Regularization,'smoothest')
  % find the smoothest possible solution
  Mfit = [0.000001*Mdes;Mreg];
  rhsfit = [0.000001*rhs;rhsreg];
else
  Mfit = [Mdes;RP*Mreg];
  rhsfit = [rhs;RP*rhsreg];
end

if prescription.Verbosity > 1
  disp(' ')
  disp(['Condition number of the regression: ',num2str(cond(Mdes))])
  disp(' ')
end

if isempty(Mineq) && isempty(Meq)
  % backslash will suffice
  coef = Mfit\rhsfit;
  lambda.eqlin=[];
  lambda.ineqlin=[];
  
  solver = 'backslash';

elseif isempty(Mineq)
  % with no inequality constraints, lse is faster than
  % is lsqlin. This also allows the use of slm when
  % the optimization toolbox is not present if there
  % are no inequality constraints.
  coef = lse(Mfit,rhsfit,full(Meq),rhseq);
  lambda.eqlin=[];
  lambda.ineqlin=[];
  
  solver = 'lse';
  
else
  % use lsqlin. first, set the options
  options = optimset('lsqlin');
  if prescription.Verbosity > 1
    options.Display = 'final';
  else
    options.Display = 'off';
  end
  % the Largescale solver will not allow general constraints,
  % either equality or inequality
  options.LargeScale = 'off';
  options.Algorithm = prescription.LsqlinAlgorithm;
  
  % and solve
  [coef,junk,junk,exitflag,junk,lambda] = ...
    lsqlin(Mfit,rhsfit,Mineq,rhsineq,Meq,rhseq,[],[],[],options); %#ok
  
  % was there a feasible solution?
  if exitflag == -2
    warning('No feasible solution was found by LSQLIN. This may reflect an inconsistent prescription set.')
    coef = nan(size(Mfit,2),1);
  end
  
  solver = 'lsqlin';
  
end

if prescription.Verbosity > 0
  disp('=========================================')
  disp(['Solver chosen as:     ',solver])
  disp('=========================================')
end

% ========================================================
% ========== chose RP to hit aim rmse ====================
% ========================================================
function delta_rmse = match_rmse(RP,Mdes,rhs,Mreg,rhsreg,Meq,rhseq,Mineq,rhsineq,aim_rmse,prescription)
% returns delta rmse for a given regularization parameter

% the regularization parameter has been logged.
% exponentiate for the solve
RP = 10^RP;

% solve the system at the given RP
coef = solve_slm_system(RP,Mdes,rhs,Mreg,rhsreg, ...
  Meq,rhseq,Mineq,rhsineq,prescription);

% compute rmse
rmse = sqrt(mean((rhs - Mdes*coef).^2));

% delta rmse error from aim
delta_rmse = abs(aim_rmse - rmse);


% ========================================================
% ======= chose RP to minimize cv prediction error =======
% ========================================================
function cvpred = min_cv(RP,Mdes,rhs,Mreg,rhsreg,Meq,rhseq,Mineq,rhsineq,prescription)
% returns delta rmse for a given regularization parameter

% the regularization parameter has been logged.
% exponentiate for the solve
RP = 10^RP;

% solve the system at the given RP, dropping out one point
% at a time from the data.
n = size(Mdes,1);
uselist = true(n,1);
cvpred = 0;
for i=1:n
  ul = uselist;
  ul(i) = false;
  coef = solve_slm_system(RP,Mdes(ul,:),rhs(ul),Mreg, ...
    rhsreg,Meq,rhseq,Mineq,rhsineq,prescription);
  pred = Mdes(i,:)*coef;
  cvpred = cvpred + (pred - rhs(i)).^2;
end


% ========================================================
% ========= choose knot placement for free knots =========
% ========================================================
function rss = free_knot_obj(intknots,x,y,prescription)
% returns residual sum of squares for fmincon to work on

% insert the interior knots from the optimizer
prescription.Knots(2:(end-1)) = intknots;
slm = slmengine(x,y,prescription);

% get predictions, form residuals, square, and sum
switch prescription.Result
  case 'slm'
    rss = sum((slmeval(x,slm,0) - y).^2);
  case 'pp'
    rss = sum((ppval(slm,x) - y).^2);
end

% ========================================================
% ========= set scaling as necessary =========
% ========================================================
function [yhat,prescription,YScale,YShift] = scaleproblem(x,y,prescription)
% chooses an appropriate scaling for the problem
% there is no need to scale x, since each knot interval is already
% implicitly normalized into [0,1]. x is passed in only for
% a few parameters, such as the integral constraint, where x would
% be needed.

if strcmp(prescription.Scaling,'on')
  % scale y so that the minimum value is 1/phi, and the maximum value phi.
  % where phi is the golden ratio, so phi = (sqrt(5) + 1)/2 = 1.6180...
  % Note that phi - 1/phi = 1, so the new range of yhat is 1. (Note that
  % this interval was carefully chosen to bring y as close to 1 as
  % possible, with an interval length of 1.)
  %
  % The transformation is:
  %   yhat = y*yscale + yshift
  phi_inverse = (sqrt(5) - 1)/2;
  
  % shift and scale are determined from the min and max of y.
  ymin = min(y);
  ymax = max(y);
  
  YScale = 1./(ymax - ymin);
  if isinf(YScale)
    % in case data was passed in that is constant, then
    % the range of y is zero. No scaling need be done then.
    YScale = 1;
  end
  
  % recover the shift once the scale factor is known.
  YShift = phi_inverse - YScale*ymin;
  
  % scale y to refect the shift and scale
  yhat = y*YScale + YShift;
  
  % finally, shift/scale each part of the prescription as necessary.
  % derivative information need only be scaled of course. And since
  % YScale will always be positive, monotonicity signs
  % and curvature signs will remain unchanged.
  
  % Those constraints that will be left untouched:
  %  C2, ConcaveDown, ConcaveUp, ConstantRegion,
  %  Decreasing, Degree, EndConditions, Envelope,
  %  Increasing, InteriorKnots, Jerk, Knots,
  %  NegativeInflection, Order, Plot, PositiveInflection,
  %  Predictions, Regularization, Result, Robust,
  %  SimplePeak, SimpleValley, Verbosity, Weights
  %
  % All other constraints will potentially be affected.
  
  % ErrorBar
  if ~isempty(prescription.ErrorBar)
    prescription.ErrorBar = prescription.ErrorBar*YScale + YShift;
  end
  
  % Integral
  if ~isempty(prescription.Integral)
    % we need to know the width of the spline domain
    support = [min(x), max(x)];
    if numel(prescription.Knots) > 1
      support(1) = min(support(1),min(prescription.Knots));
      support(2) = max(support(2),max(prescription.Knots));
    end
    
    prescription.Integral = prescription.Integral* ...
      YScale + diff(support)*YShift;
  end

  % LeftMaxSlope
  if ~isempty(prescription.LeftMaxSlope)
    prescription.LeftMaxSlope = prescription.LeftMaxSlope*YScale;
  end
  
  % LeftMaxValue
  if ~isempty(prescription.LeftMaxValue)
    prescription.LeftMaxValue = prescription.LeftMaxValue*YScale + YShift;
  end
  
  % LeftMinSlope
  if ~isempty(prescription.LeftMinSlope)
    prescription.LeftMinSlope = prescription.LeftMinSlope*YScale;
  end
  
  % LeftMinValue
  if ~isempty(prescription.LeftMinValue)
    prescription.LeftMinValue = prescription.LeftMinValue*YScale + YShift;
  end
  
  % LeftSlope
  if ~isempty(prescription.LeftSlope)
    prescription.LeftSlope = prescription.LeftSlope*YScale;
  end
  
  % LeftValue
  if ~isempty(prescription.LeftValue)
    prescription.LeftValue = prescription.LeftValue*YScale + YShift;
  end
  
  % LinearRegion
  % nothing to be done for this one
  
  % MaxSlope
  if ~isempty(prescription.MaxSlope)
    prescription.MaxSlope = prescription.MaxSlope*YScale;
  end

  % MaxValue
  if ~isempty(prescription.MaxValue)
    prescription.MaxValue = prescription.MaxValue*YScale + YShift;
  end
  
  % MaxFpp
  if ~isempty(prescription.MaxFpp)
    prescription.MaxFpp = prescription.MaxFpp*YScale;
  end
  
  % MinSlope
  if ~isempty(prescription.MinSlope)
    prescription.MinSlope = prescription.MinSlope*YScale;
  end

  % MinValue
  if ~isempty(prescription.MinValue)
    prescription.MinValue = prescription.MinValue*YScale + YShift;
  end
  
  % MinFpp
  if ~isempty(prescription.MinFpp)
    prescription.MinFpp = prescription.MinFpp*YScale;
  end
  
  % RightMaxSlope
  if ~isempty(prescription.RightMaxSlope)
    prescription.RightMaxSlope = prescription.RightMaxSlope*YScale;
  end
  
  % RightMaxValue
  if ~isempty(prescription.RightMaxValue)
    prescription.RightMaxValue = prescription.RightMaxValue*YScale + YShift;
  end
  
  % RightMinSlope
  if ~isempty(prescription.RightMinSlope)
    prescription.RightMinSlope = prescription.RightMinSlope*YScale;
  end

  % RightMinValue
  if ~isempty(prescription.RightMinValue)
    prescription.RightMinValue = prescription.RightMinValue*YScale + YShift;
  end
  
  % RightSlope
  if ~isempty(prescription.RightSlope)
    prescription.RightSlope = prescription.RightSlope*YScale;
  end
  
  % RightValue
  if ~isempty(prescription.RightValue)
    prescription.RightValue = prescription.RightValue*YScale + YShift;
  end
  
  % Scaling ...
  % once y has been shifted and scaled, we need not do it again
  prescription.Scaling = 'off';
  
  % SumResiduals
  if ~isempty(prescription.SumResiduals)
    prescription.SumResiduals = prescription.SumResiduals*YScale;
  end
  
  % XY
  if ~isempty(prescription.XY)
    prescription.XY(:,2) = prescription.XY(:,2)*YScale + YShift;
  end
  
  % XYP
  if ~isempty(prescription.XYP)
    prescription.XYP(:,2) = prescription.XYP(:,2)*YScale;
  end
  
  % XYPP
  if ~isempty(prescription.XYPP)
    prescription.XYPP(:,2) = prescription.XYPP(:,2)*YScale;
  end
  
  % XYPPP
  if ~isempty(prescription.XYPPP)
    prescription.XYPPP(:,2) = prescription.XYPPP(:,2)*YScale;
  end
  
else
  % no scaling is done
  YShift = 0;
  YScale = 1;
  % so we will leave y and the prescription unchanged
  yhat = y;
end

prescription.YScale = YScale;
prescription.YShift = YShift;

% ========================================================
% ========= chooseknots =========
% ========================================================
function [knots,nk] = chooseknots(K,n,x)
% choose every K'th data point as a knot

if K > 0
  % just given a number of knots, equally spaced
  nk = K;
  knots = linspace(min(x),max(x),nk)';
  knots(end)=max(x); % just to make sure
else
  % we need every abs(K)'th knot
  
  if (K == 0) || (K >= n)
    error('SLMENGINE:knots','Every K''th data point was specified, but K was too large or zero')
  end
  
  Kind = 1:abs(K):n;
  % in case K was not an integer, or did not
  % go to the very end.
  if (n - Kind(end)) < (K/2)
    % expand the last knot interval
    Kind(end) = n;
  else
    % append one extra knot for the last data point
    Kind(end + 1) = n;
  end
  % rounding makes them into indices
  Kind = round(Kind);
  % in case of replicate data points
  x = sort(x);
  knots = unique(x(Kind));
  nk = numel(knots);
end



% ========================================================
% ========= model statistics =========
% ========================================================
function slmstats = modelstatistics(slmstats,y,coef,YScale,Mdes)

% generate model statistics, stuffing them into slmstats

% residuals, as yhat - y
resids = (Mdes*coef - y)./YScale;

% RMSE: Root Mean Squared Error
slmstats.RMSE = sqrt(mean(resids.^2))./YScale;

% R-squared
slmstats.R2 = 1 - sum(resids.^2)./sum(((y - mean(y))./YScale).^2);

% adjusted R^2
ndata = numel(y);
slmstats.R2Adj = 1 - (1-slmstats.R2)*(ndata - 1)./(ndata - slmstats.NetDoF);

% range of the errors, min to max, as yhat - y
slmstats.ErrorRange = [min(resids),max(resids)];

% compute the 25% and 75% points (quartiles) of the residuals
% (This is consistent with prctile, from the stats TB.)
resids = sort(resids.');
ind = 0.5 + ndata*[0.25 0.75];
f = ind - floor(ind);
ind = min(ndata - 1,max(1,floor(ind)));
slmstats.Quartiles = resids(ind).*(1-f) + resids(ind+1).*f;


% ========================================================
% ========== check for inappropriate properties ==========
% ========================================================
function property_check(prescription,model_degree)
% issues warning messages for inappropriate properties
% for the given model degree.

if strncmpi(prescription.Robust,'on',2)
  warning('SLMENGINE:RobustFitting', ...
    'Sorry, but the robust fitting option is not yet operable.')
  prescription.Robust = 'off';
end

switch model_degree
  case {3 'cubic'}
    % no properties flagged for cubic models
    
  case {1 'linear'}
    % only (some) curvature properties are flagged
    % for the linear model
    
    if ~isempty(prescription.Jerk)
      warning('SLMENGINE:ignoredconstraint', ...
        'Property ignored for a linear model: Jerk')
    end
    
    if ~isempty(prescription.MinFpp)
      warning('SLMENGINE:ignoredconstraint', ...
        'Property ignored for a linear model: MinFpp')
    end
    
    if ~isempty(prescription.MaxFpp)
      warning('SLMENGINE:ignoredconstraint', ...
        'Property ignored for a linear model: MaxFpp')
    end
    
    if ~isempty(prescription.SegmentLinear)
      warning('SLMENGINE:ignoredconstraint', ...
        'Property ignored for a linear model: SegmentLinear')
      prescription.SegmentLinear = [];
    end
    
    if ~isempty(prescription.XYPP)
      warning('SLMENGINE:ignoredconstraint', ...
        'Property ignored for a linear model: XYPP')
    end
    
    if ~isempty(prescription.XYPPP)
      warning('SLMENGINE:ignoredconstraint', ...
        'Property ignored for a linear model: XYPPP')
    end
    
    if strncmpi(prescription.Extrapolation,'cubic',5)
      warning('SLMENGINE:extrapolation', ...
        'Linear model will not employ cubic extrapolation')
      prescription.Extrapolation = 'linear';
    end
    
    
  case {0 'constant'}
    % curvature properties are flagged for
    % the constant model, as well as most slope
    % related properties
    
    if strncmpi(prescription.Extrapolation,'cubic',5)
      warning('SLMENGINE:extrapolation', ...
        'Constant model will not employ cubic extrapolation')
      prescription.Extrapolation = 'constant';
    end
    
    if strncmpi(prescription.Extrapolation,'linear',6)
      warning('SLMENGINE:extrapolation', ...
        'Constant model will only employ constant extrapolation')
      prescription.Extrapolation = 'constant';
    end
    
    if ~isempty(prescription.Jerk)
      warning('SLMENGINE:ignoredconstraint', ...
        'Property ignored for a constant model: Jerk')
    end
    
    if ~isempty(prescription.MinFpp)
      warning('SLMENGINE:ignoredconstraint', ...
        'Property ignored for a constant model: MinFpp')
    end
    
    if ~isempty(prescription.MaxFpp)
      warning('SLMENGINE:ignoredconstraint', ...
        'Property ignored for a constant model: MaxFpp')
    end
    
    if ~isempty(prescription.LeftMaxSlope)
      warning('SLMENGINE:ignoredconstraint', ...
        'Property ignored for a constant model: LeftMaxSlope')
    end
    
    if ~isempty(prescription.LeftMinSlope)
      warning('SLMENGINE:ignoredconstraint', ...
        'Property ignored for a constant model: LeftMinSlope')
    end
    
    if ~isempty(prescription.LeftSlope)
      warning('SLMENGINE:ignoredconstraint', ...
        'Property ignored for a constant model: LeftSlope')
    end
    
    if ~isempty(prescription.LinearRegion)
      warning('SLMENGINE:ignoredconstraint', ...
        'Property ignored for a constant model: LinearRegion')
    end
    
    if ~isempty(prescription.MaxSlope)
      warning('SLMENGINE:ignoredconstraint', ...
        'Property ignored for a constant model: MaxSlope')
    end
    
    if ~isempty(prescription.MinSlope)
      warning('SLMENGINE:ignoredconstraint', ...
        'Property ignored for a constant model: MinSlope')
    end
    
    if ~isempty(prescription.NegativeInflection)
      warning('SLMENGINE:ignoredconstraint', ...
        'Property ignored for a constant model: NegativeInflection')
    end
    
    if ~isempty(prescription.PositiveInflection)
      warning('SLMENGINE:ignoredconstraint', ...
        'Property ignored for a constant model: PositiveInflection')
    end

    if ~isempty(prescription.RightMaxSlope)
      warning('SLMENGINE:ignoredconstraint', ...
        'Property ignored for a constant model: RightMaxSlope')
    end
    
    if ~isempty(prescription.RightMinSlope)
      warning('SLMENGINE:ignoredconstraint', ...
        'Property ignored for a constant model: RightMinSlope')
    end
    
    if ~isempty(prescription.RightSlope)
      warning('SLMENGINE:ignoredconstraint', ...
        'Property ignored for a constant model: RightSlope')
    end
    
    if ~isempty(prescription.SegmentConstant)
      warning('SLMENGINE:ignoredconstraint', ...
        'Property ignored for a constant model: SegmentConstant')
      prescription.SegmentConstant = [];
    end
    
    if ~isempty(prescription.SegmentLinear)
      warning('SLMENGINE:ignoredconstraint', ...
        'Property ignored for a constant model: SegmentLinear')
      prescription.SegmentLinear = [];
    end
    
    if ~isempty(prescription.XYP)
      warning('SLMENGINE:ignoredconstraint', ...
        'Property ignored for a constant model: XYP')
    end
    
    if ~isempty(prescription.XYPP)
      warning('SLMENGINE:ignoredconstraint', ...
        'Property ignored for a constant model: XYPP')
    end
    
    if ~isempty(prescription.XYPPP)
      warning('SLMENGINE:ignoredconstraint', ...
        'Property ignored for a constant model: XYPPP')
    end
    
end



